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ABSTRACT 


A  model  for  a  dam  is  considered  wherein  the  net  input  rate  (input 
minus  output  rate)  follows  a  finite  Markov  chain  in  continuous  time, 

Xt,  and  the  dam  contents  process,  C^.,  is  the  integral  of  the  Markov 
chain.  The  dam  is  then  modelled  with  the  bivariate  Markov  process 
(Xt,Ct),  of  which  three  variations  are  considered.  These  are  the 
doubly-infinite  dam  with  no  top  or  bottom,  the  semi-infinite  dam  with 
only  one  boundary,  and  the  finite  dam  with  both  a  top  and  a  bottom. 

Some  of  the  analysis  is  performed  under  the  most  general  situation  in 
which  X^  is  defined  on  m  states  and  has  an  arbitrary  generator,  while 
other  analysis  is  performed  under  the  restricted  case  when  m  =  2. 

For  the  doubly-infinite  dam,  the  first  and  second  moment  functions 
and  the  maximum  and  minimum  variables  are  studied.  The  expected  range 
function  is  explicitly  derived  in  a  special  two-state  case.  Also  in 
the  two-state  case,  weak  convergence  to  the  Wiener  process  is  established 
in  D[0/»),  from  which  the  asymptotic  distribution  of  the  range  is  obtained. 

For  the  semi-infinite  and  finite  dams,  techniques  of  invariance  used 
in  the  physical  sciences  are  introduced  to  study  first  passage  times. 

These  techniques  are  used  to  derive  the  distribution  and  moments  of  the 
wet  period  of  the  dam  in  special  cases,  and  the  limiting  probabilities 
of  emptiness  and  overflow. 
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CHAPTER  I 


INTRODUCTION 


1 .1  Prel iminaries 

The  field  of  storage  theory  has  been  of  long-standing  interest 
in  the  engineering  community,  and  has  in  the  last  few  decades  become 
one  of  the  most  energetically  pursued  fields  in  applied  probability. 
Although  the  origin  of  storage  theory  lies  in  the  now  classical  study 
of  sizing  of  water  reservoirs  by  hydrologists,  the  recent  developments 
have  transcended  these  original  applications.  The  interest  of 
probabilists  and  statisticians  in  the  interesting  mathematical  and 
statistical  problems  which  have  arisen  from  a  study  of  storage  problems 
have  in  fact  led  to  a  new  branch  in  hydrology  known  as  stochastic 
hydrology.  Many  of  the  techniques  developed  by  the  stochastic 
hydrologists  have  in  turn  become  of  independent  interest,  so  that  at 
present  this  theory  has  developed  on  its  own  merits  as  a  mathematical 
construct,  rather  than  being  strictly  tied  to  the  classical  study  of 
reservoir  sizing.  Many  natural  linkages  have  also  been  established 
between  storage  theory  and  other  areas  in  stochastic  processes,  among 
them  queuing  theory  and  the  theory  of  stationary  processes. 

As  a  result  of  the  extensiveness  of  the  field,  it  is  necessary 
for  contemporary  investigators  to  narrow  substantially  the  field  of 
inquiry  to  specific  types  of  models  and  specific  aspects  of  those 
models.  This  investigation  studies  certain  problems  associated  with 
a  continuous-time  storage  model  which  has  a  Markovian  structure. 

8efore  specifically  addressing  the  problems  studied  in  this 
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investigation,  however,  we  will  give  a,  necessarily  brief,  overview 
of  the  history  and  categorizations  of  storage  theory. 

1 .2  History  and  Literature  Review 

The  classic  object  of  storage  theory  is  the  study  of  a  water 
reservoir,  or  dam,  which  is  fed  by  some  source  of  water  such  as  a 
river.  The  size  of  the  dam  to  be  constructed  is,  of  course,  dependent 
on  the  amount  of  water  which  can  be  expected  from  the  supply  source. 

The  size  of  the  dam  should  be  such  that  the  probability  of  occurrence 
of  either  overflow  or  emptiness  is  minimized. 

The  classical  analysis  of  this  problem  specified  certain 
deterministic  functions  for  both  the  amount  of  water  input  and  the 
amount  of  water  drawn  off,  and  from  this  the  optimal  size  of  the 
reservoir  was  determined.  This  type  of  analysis  dates  at  least  as 
far  back  as  1883,  when  it  was  treated  by  W.  Rippl. 

The  stochastic  nature  of  the  problem  was  first  addressed  by  A. 

Hazen  (1914)  and  later  by  C.  E.  Sudler  (1927).  These  two  individuals 
considered  the  data  representing  river  runoffs  as  having  an  associated 
uncertainty.  Although  their  methods  of  analysis  were  crude  by  today's 
standards,  it  is  both  to  their  credit  and  an  indication  of  the 
interesting  nature  of  the  problems  in  storage  theory  that  they 
developed  the  techniques  of  using  "probability  paper"  and  data 
simulation,  respectively,  in  their  studies. 

H.  E.  Hurst  [26]*  in  1951  studied  the  problem  of  determining 
the  reservoir  storage  required  on  a  given  stream  to  guarantee  a  given  • 

draft  by  considering  the  cumulative  sums  of  the  departure  of  the 
♦Numbers  in  brackets  refer  to  bibliography. 
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annual  totals  from  the  mean  annual  total  discharge.  These  cumulative 

departures  can  be  thought  of  as  the  contents  of  a  hypothetical  dam 

with  no  top  or  bottom;  or  the  unrestricted  contents,  measured  from  s.ne 

point  representing  the  mean  contents. 

Under  this  model,  the  level  of  the  dam  in  year  n  ,  Cn  ,  is 

the  result  of  n  net  yearly  inputs  (input  minus  output)  X^  X2,  ... 

n 

Xn  .  Thus  we  have  Cn  s[  X^  ,  so  that  the  contents  become  the  well- 
studied  partial  sums  process.  The  functional  studied  by  Hurst  for 
this  model  is  the  range  up  to  time  n,  R  ,  given  by 


p  _  max  C-  min  C • 

n  1  <  i  <  n  1  <  i  <  n 


Hurst  and  subsequent  investigators  used  Rn  for  initial  sizing 
estimates  by  concluding  that  the  range  up  to  time  n  gives  an 
indication  of  the  size  of  a  dam  which  would  have  been  required  to 
contain  that  amount  of  water  without  either  overflowing  or  becoming 
empty.  Hurst  also  considered  certain  other  variations  on  the  range 
as  defined  above,  such  as  an  adjusted  range  obtained  by  subtracting 
off  the  average  up  to  time  n  .  The  basic  idea  for  the  sizing  of  the 
reservoir,  however,  is  the  same. 

If  the  interval  size  in  this  discrete-time  model  is  sufficiently 
large  (yearly,  for  example),  then  it  may  be  reasonable  to  suppose  that 
the  summands  X..  are  stochastically  independent.  This  was  the  first 
model  studied,  and  in  this  case  the  unrestricted  contents  process 
becomes  a  random  walk. 

The  significance  of  the  Hurst  paper  lies  in  the  extremely  long 
records  of  data  which  he  compiled,  some  extending  up  to  2000  years,  as 
was  the  case  for  the  records  on  the  Nile  river.  From  these  records  he 
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drew  a  conclusion  on  the  rate  of  growth  of  the  range  which  contra¬ 
dicted  the  random  walk  model  described  above.  For  a  random  walk 
model  it  is  well  known  that  the  expected  value  of  the  range  grows  as 
n^‘5  .  However,  Hurst  found  in  his  analysis  of  the  data  that  all  of 
these  quantities  showed  a  surprisingly  similar  exponent  of  0.69  to 
0.80  with  a  mean  of  .72.  This  is  indeed  much  too  large  to  be  explained 
by  this  model  and  has  been  the  object  of  much  subsequent  study. 

This  anomaly,  now  known  as  the  Hurst  phenomenon,  drew  the 
attention  of  W.  Feller  [21]  who  in  1951  derived  the  asymptotic 
distribution  and  moments  of  the  range  for  the  iid  case  by  appealing 
to  approximations  by  Brownian  motion.  Feller  mentioned  that  the 
Hurst  phenomenon  might  be  explained  by  assuming  that  the  summands 
are  not  independent,  although  as  noted  later  by  P.  Moran  [33],  the 
dependence  would  have  to  be  strange  indeed,  for  any  reasonable  model 
will  in  fact  have  an  asymptotic  growth  of  0.5.  Moran  then  commented 
that  the  Hurst  behavior  may  be  pre-asymptotic  in  nature,  and  that  the 
records  which  Hurst  studied  were  not  long  enough  to  reach  their 
asymptotic  values. 

Subsequently  Moran  in  1964  [34]  obtained  the  mean  range  when  the 
inputs  were  iid  but  with  a  symmetric  stable  distribution  with 
parameter  y  and  found  that  the  mean  range  varied  as  n^Y. 

This,  then,  represents  the  two  main  theories  advanced  to  the 
present  time  to  explain  the  Hurst  phenomenon;  non-iid  pre-asymptotic 
behavior,  and  iid  heavy-tailed  net  inputs  as  represented  by  the  stable 
inputs. 

The  heavy  tails  explanation  is  not  as  appealing  to  many  applied 
hydrologists  because  of  conceptual  difficulties  involved  with  an  input 
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which  has  an  infinite  variance.  However,  it  is  the  opinion  of  this 
investigator  that  the  pre-asymptotic  theory  is  burdened  by  the  fact 
that  none  of  the  Hurst  data  seemed  asymDtotic  tc  n  *  .  In  fact,  a1!  of 

the  series  which  he  examined  showed  remarkably  consistent  large 
growth  rates  for  as  far  back  as  he  could  find  data. 

Although  we  will  not  concern  ourselves  in  great  detail  with  the 
Hurst  behavior,  we  will  demonstrate  that  the  continuous- time  model 
which  we  investigate  does  in  fact  exhibit  the  Hurst  behavior  pre- 
asymptotical ly  for  moderate  values  of  time. 

The  more  modern  studies  in  storage  theory  began  with  Moran  in  a 
series  of  papers  from  1954  -  1957  in  which  he  studied  dams  with  a 
variety  of  inputs  and  operational  policies.  In  his  analysis,  Moran 
considered  the  finite  dam  with  both  a  top  and  bottom  and  performed 
the  analysis  in  both  discrete  time  and  continuous  time. 

In  the  continuous  time  approach,  Moran  assumed  that  the  input 
process  was  an  additive  homogeneous  process;  that  is,  a  process  with 
stationary  independent  increments.  This  means  that  the  input  incre¬ 
ments,  say  C{t2)  -  C ( 1 1 )  in  non-overlapping  intervals  (t1,t2)  are 
independent  and  have  a  distribution  which  depends  only  on  t2  -  t1  . 

This  is  a  straight  analog  of  the  situation  in  discrete  time  when  the 
successive  inputs  are  considered  to  be  independent.  Although  Moran’s 
continuous  time  analysis  is  a  more  convenient  model  for  a  real  dam  in 
that  one  is  no  longer  limited  to  studying  the  inputs  only  at  certain 
time  epochs,  it  nevertheless  has  the  fault  that  as  the  time  intervals 
become  shorter,  the  assumption  of  independence  in  successive  intervals 
becomes  less  realistic.  Among  other  contributors  to  the  model  is 
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J.  Gani  [23],  who  considered  the  input  process  to  be  Poisson  with  a 
unit  release  rate. 

E.  Lloyd  [30]  in  1963  was  the  first  to  consider  a  dam,  in 
discrete  time,  where  the  successive  inputs  were  dependent.  He 
supposed  that  the  input  Xn  in  successive  time  intervals  (n,  n+1 ) 
followed  a  Markov  chain  and  the  release  rate  in  each  interval  was  a 
constant,  say  r  .  The  dam  contents  at  time  n  are  then  given  by 
Cn  =  min (a,  max{0,  Cn_x  +  Xn-r)) 
where  a  is  the  capacity  of  the  dam  and  0  <_  XQ  <_  a  . 

This  model,  which  has  since  become  known  as  the  Lloyd  dam, 
received  considerable  attention  in  the  discrete  time  case  in  both  this 
version  and  the  semi-infinite  topless  version  where  a  *  °°  .  Notable 
contributions  to  this  study  include  Ali  Khan  and  Gani  [1],  who  studied 
the  time  to  first  emptiness  for  the  semi -infinite  dam,  and  Ali  Khan 
[2],  who  considered  the  finite  case.  In  all  of  these  studies  the 
Markov  chain  has  a  finite  state  space.  Brockwell  and  Garti  [11]  con¬ 
sidered  the  time  to  first  emptiness  for  the  case  in  which  the  Markov 
chain  has  the  non-negative  integers  as  state  space. 

As  far  as  range  analysis  for  the  unrestricted  contents  with 
dependence  is  concerned,  F.  Gomide  [24]  treated  the  case  of  Markovian 
inputs  in  discrete  time,  and  B.  Troutman  [41]  studied  limiting  dis¬ 
tributions  for  the  discrete  time  process  with  Markovian  and  certain 
stationary  inputs  using  a  weak  convergence  approach. 

If  the  continuous  time  model  of  a  storage  system  is  a  better 
approximation  than  the  corresponding  discrete  time  model,  and  the 
Lloyd  model  is  a  better  approximation  than  the  independent  input  model, 
then  a  continuous  time  version  of  the  Lloyd  dam  offers  a  much  closer 
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correspondence  between  theory  and  a  realistic  storage  system. 
Relatively  few  studies  have  been  done  on  this  model,  which  is  the 
subject  of  this  investigation.  Most  notable  of  the  published  result; 
have  been  those  by  McNeil  [32]  and  Brockwell  [12].  McNeil  studied  a 
dam  in  continuous  time  in  which  the  input  process  follows  a  two-state 
Markov  chain  with  one  of  the  states  being  zero.  When  the  input 
process  is  in  the  zero  state,  the  dam  is  being  drained,  and  he 
supposes  that  there  is  a  general  measurable  function  of  the  dam  level, 
say  g(x),  which  represents  the  demand  rate.  Although  he  sets  up  the 
problem  in  this  generality,  he  is  able  to  obtain  explicit  results 
only  for  the  cases  when  g(x)  is  constant  and  g(x)  is  exponential. 
McNeil  is  able  to  derive  the  limiting  distribution  of  the  contents  and 
first  passage  times. 

Brockwell  considers  the  case  in  which  the  net  input  rate  follows 
a  general  Markov  chain  with  a  finite  state  space,  and  by  setting  up 
Kolmogorov-type  equations  is  able  to  derive  the  limiting  distribution 
of  the  contents  and  first  times  to  emptiness  and  overflow. 

1 .3  Objectives  and  General  Approach 

The  model  considered  in  this  investigation  is  the  one  formulated 
by  Brockwell  where  the  net  input,  say  Xt  ,  is  a  Markov  chain  in 
continuous  time  defined  on  a  finite  state  space  {u1,...,u|T)}  .  In 
analogy  with  discrete-time  models,  the  dam  content  at  time  t  ,  Ct  , 
is  given  by  the  integral  of  the  net  input  process  up  to  time  t  . 

Hence  (for  a  doubly  infinite  dam) 
t 

Ct  «  /  xudu,  t  _>  0  . 
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This  formulation  is  the  continous-time  analog  of  the  unrestricted 

contents  process  described  earlier.  To  contrast  it  with  the  next  model 

which  we  consider,  we  will  refer  to  the  bivariate  process  (X t )  as 

the  doubly-infinite  dam.  When  we  refer  to  the  marginal  process  C  , 

the  term  "unrestricted  contents"  will  also  be  used. 

If  we  now  refine  this  model  by  assuming  that  the  contents  must  be 

non-negative  then  £  .  j*  x,dllj  t  „  ^,ere 

t  0  u 

X*  =  0  if  (J  -  o  and  Xy  <  0 

u  Xu  otherwise 

and  we  will  refer  to  the  bivariate  process  (Xt,cJ)  in  this  context  as 
the  semi -infinite  dam.  The  terminology  topless  dam  has  also  been  used 
for  models  of  this  type.  A  symmetric  variation  of  this  is  to  restrict 
the  top  but  not  the  bottom,  so  that,  if  a  is  the  highest  level  which 
the  contents  can  attain,  then  3  =  ft  ^  t  Q  §  where 

^  O  ^ 

u  B 

X*  _  0  if  Cu  =  a  and  Xu  >  0 

u  Xu  otherwise. 

We  will  refer  to  this  variation  of  the  semi-infinite  dam  as  the 
bottomless  dam. 

The  third  model  which  we  consider  is  one  for  which  the  contents 
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Now,  in  simple  terms  a  Markov  process  is  one  in  which  the 
distribution  of  the  future,  given  the  present,  is  independent  of  the 
past.  The  Ct  process  does  not  meet  the  requirement  and  is  not,  in 
fact,  Markov.  This  is  due  to  the  fact  that  the  present  only  conveys 
information  on  the  level  of  the  dam  contents,  and  not  on  the  rate  at 
which  the  contents  are  changing.  Obviously,  adding  this  information 
will  affect  the  future  distribution  of  the  contents.  The  bivariate 
process  (Xt,  Ct)  ,  however,  does  contain  the  information  on  the  rate 
of  change  of  the  contents  and  is  a  bivariate  Markov  process.  The 
state  space  of  such  a  process  is  discontinuous  i  t  nature,  since  it 
evolves  on  a  set  of  disjoint  lines  as  shown  in  Figure  1.1  for  the 
finite  dam.  The  solid  dots  at  the  right  boundaries  for  the  positive 
rates  and  the  left  boundaries  for  the  negative  rates  represent  the 
fact  that  when  the  process  hits  these  points  it  remains  there  for  a 
random  time.  A  solid  dot  on  the  right  boundary,  for  example, 
represents  an  overflow  condition  while  a  solid  dot  on  the  left 
boundary  represents  an  emptiness  condition.  Markov  processes  with 
discontinuous  state  spaces  of  this  type  were  studied  by  Moyal  [35]. 

In  Chapter  II  we  present  a  brief  review  of  background  material 
which  is  prerequisite  to  the  development  of  later  chapters.  Although 
none  of  the  material  is  new,  the  results  which  we  quote  here  are 
scattered  through  various  sources  and  so  we  include  it  here  in  a 
unified  manner  as  an  aid  to  the  reader. 

In  Chapter  III  we  study  the  first  and  second  moment  functions  of 
the  marginal  processes  for  the  doubly  infinite  dam.  These  functions 
are  useful  for  fitting  the  model  and  as  an  aid  in  the  initial  estimates 
of  parameters.  In  this  chapter  we  prove  a  generalization  of  a  result 
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reported  by  McNeil  [32]  that  the  autocorrelation  function  for  the 
marginal  contents  process  does  not  depend  on  the  states  of  the 
Markov  chain  in  the  two-state  case. 

In  Chapter  IV  we  begin  with  the  analysis  of  the  range  for  the 
doubly  infinite  dam.  The  relevance  of  this  to  dam  sizing  studies 
has  been  discussed  earlier.  Although  the  range  variable  has  been 
studied  extensively  in  discrete  time,  this  is  the  first  analysis  in 
continuous  time.  We  first  study  the  joint  distribution  of  the 
maximum  and  minimum  variables  by  finding  its  Laplace  transform.  We 
do  this  for  the  most  general  situations  in  which  the  Markov  chain  is 
defined  on  m  states  with  an  arbitrary  generator.  We  call  this 
the  General  Case.  An  inversion  of  this  transform  would  provide  the 
joint  distribution  from  which  the  distribution  of  the  range  could,  in 
principle,  be  obtained.  The  inverse,  however,  is  not  obtainable  in  a 
simple  analytic  form  and  although  a  numerical  procedure  could  be 
used,  we  choose  to  restrict  ourselves  to  an  analytic  rather  than 
numerical  analysis  of  this  problem.  By  restricting  ourselves  to 
special  cases,  we  are  able  to  obtain  exact  expressions  for  the  expected 
range  for  all  t  ,  and  an  asymptotic  distribution  for  the  range.  The 
special  cases  which  we  consider  are  those  in  which  the  Markov  chain 
Is  defined  cn  two  states,  which  we  call  the  Two-State  Case  and  the 
subcase  for  which  the  two  states  ui  and  u2  ,  and  the  two  holding 
time  parameters,  A  and  p  ,  have  the  relation  m  =  -u2  and  A  *  p  . 
We  call  this  the  Symmetric  Case. 

We  are  able  to  derive  an  explicit  solution  for  the  expected  range 
function,  ERt  ,  in  the  symmetric  case.  We  do  this  by  first  obtaining 
the  marginal  distributions  of  the  maximum  and  minimum  for  the  two-state 


. •  .1.  .. -yw*  .  * 
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case.  The  expected  range  can  be  expressed,  using  the  linearity  of 
expectation,  as  the  difference  in  expected  value  of  the  maximum  and 
minimum.  We  obtain  the  expected  range  function  in  its  Laplace 
transformed  version.  Fortunately,  the  Laplace  transform  can  be  expanded 
in  a  neighborhood  of  infinity  and  the  resulting  series  can  be  inverted 
term-by-term  as  a  confluent  hyper geometric  series.  The  properties 
of  the  confluent  hypergeometric  series  have  been  thoroughly  studied 
(see,  e.g.,  Buchholz  [15]).  Tabled  values  are  available  in  Jahnke, 

Emde,  and  Losch  [27]. 

For  illustrative  purposes,  we  chose  a  particular  subcase  and 
performed  a  precise  numerical  calculation  of  ERt  which  we  include  as 
a  graph  for  0  <_  t  <_  10  in  figure  4.1.  For  numerical  calculation 
purposes,  Kummer's  first  formula  (Buchholz)  was  quite  useful,  for  it 
allowed  calculation  in  positive  terms  rather  than  in  terms  of  a  slowly 
converging  alternating  series. 

Examination  of  the  graph  of  ER^  shows  remarkably  close  behavior 
to  that  reported  by  Hurst  for  moderate  values  of  t  .  Therefore  we 
include  in  this  chapter  a  short  note  on  the  Hurst  phenomenon.  Another 
aspect  which  is  evident  from  the  graph  is  the  extremely  rapid  con¬ 
vergence  of  ERt  to  its  asymptotic  value  of  /\  This  indicates 

that  for  moderate  values  of  t  ,  as  small  as  5,  very  good  approxima¬ 
tions  can  be  obtained  from  asymptotic  results.  In  conclusion,  we 
discuss  the  order  of  convergence  by  performing  an  asymptotic  expansion 
of  the  Laplace  transform  of  ERt  in  a  neighborhood  of  zero,  as 
discussed  in  Doetsch  [19]. 

In  Chapter  V  we  discuss  what  we  call,  in  general,  invariance 
methods  as  applied  to  the  analysis  of  the  semi-infinite  and  the  finite 
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dams.  The  invariance  techniques  to  which  we  refer  are  those  developed 
by  V.  A.  Ambarzumian  [3]  and  Chandrasekhar  [16]  in  the  1 940 * s  for  the 
study  of  radiative  transfer  in  stellar  atmospheres,  and  subsequently 
refined  by  Bellman  et  at.  [4-7]  in  the  late  195Q‘s,  and  early  1 960 ‘ s 
in  the  field  of  Neutron  Transport  Theory.  Invariance  techniques  are 
very  closely  related  to  the  idea  of  regeneration  in  probability  theory. 
In  this  chapter  we  present  a  detailed  discussion  of  the  principle  of 
invariance  as  applicable  to  our  model. 

The  principle  of  invariance  presented  the  astrophysicists  with 
a  powerful  tool  for  the  solution  of  certain  physical  problems  which 
had  been  at  best  laboriously  solved  using  more  classical  techniques. 

The  same  held  true  in  the  applications  to  Neutron  Transport  Theory. 

In  studying  the  problems  of  first  passage  times  in  the  storage 
model  considered  herein,  a  strong  relation  between  all  of  these 
problems  becomes  evident.  This  is  particularly  true  in  the  study  of 
the  wet  period,  or  the  sojourn  time  of  the  process  from  the  time  that 
it  leaves  the  zero  state  until  the  time  that  it  returns  to  the  zero 
state.  The  classical  formulation  for  the  distribution  of  this  variable 
involves  setting  up  Kolmogorov  equations,  which  lead  to  a  two-point 
boundary  value  problem  in  the  finite  dam.  The  invariance  techniques, 
however,  enable  us  to  solve  directly  for  the  transform  of  this 
distribution,  leading  to  a  system  of  algebraic  equations  in  the  case  of 
the  semi-infinite  dam,  and  to  an  initial  value  problem  in  the  case  of 
the  finite  dam.  In  the  symmetric  case,  we  are  able  to  invert  the  trans¬ 
form  and  thus  obtain  the  distribution  of  the  wet  period.  More 
generally,  we  are  able  to  solve  for  the  expected  value  of  the  wet  ; 
period  in  the  two-state  case,  which  in  turn  leads  us  to  finding 
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necessary  and  sufficient  conditions  for  recurrence  for  the  semi- 
infinite  dam. 

We  finish  this  chapter  by  using  a  renewal  argument  to  calculate 
the  limiting  probability  of  emptiness  for  the  topless  and  finite  dams, 
and  the  limiting  probability  of  overflow  for  the  bottomless  dam. 

In  Chapter  VI  we  establish  a  Functional  Central  Limit  Theorem 
for  the  doubly  infinite  dam  on  two  states.  We  prove  that  the  process 
C(nt)/*/n  converges  weakly  in  the  zero  drift  case  to  the  Wiener  process 
as  n  goes  to  infinity,  that  is,  C(n* )//rT  =->  W(B* )  on  D[0,°»)  where 
6  is  a  function  of  the  four  parameters.  We  therefore  have  the  approxi¬ 
mation  P[C(n-)///n  eB]  =  P[W(0-)eB]  ,  for  BeB  ,  the  Borel  field  induced 
by  the  open  sets  relative  to  the  D[0,“>)  metric.  This  is  a  sub¬ 
stantial  improvement  on  the  results  of  Fukushima  and  Hitsuda  [22]  and 
Pinsky  [38],  who  were  only  able  to  show  convergence  of  the  marginal 
distributions. 

We  establish  weak  convergence  by  the  method  of  establishing  an 
embedded  partial  sum  process.  The  difficulty  here  is  that  a  direct 
appeal  to  Donsker's  Theorem  is  not  possible  because  the  sum  contains  a 
renewal  counting  function  as  an  upper  index.  Fortunately,  we  are  able 
to  overcome  this  difficulty  by  a  technique  similar  to  that  used  by 
Resnick  and  Durrett  [20]  who  consider  weak  convergence  of  sums  with 
random  indices. 

As  an  application  of  the  weak  convergence  which  we  establish,  we 
use  the  continuous  mapping  theorem  to  establish  the  asymptotic 
distribution  of  the  range,  since  the  corresponding  result  for  the 
Wiener  process  was  established  by  Feller  [21]. 
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Chapter  VII  presents  a  summary  of  the  results  and  recommendations 
for  the  further  study  of  these  problems. 


CHAPTER  II 


BACKGROUND  MATERIAL 


2.1  Markov  Processes 

In  this  chapter  we  present  a' very  brief  summary  of  definitions 
and  results  which  will  be  needed  in  subsequent  chapters.  Proofs  are 
not  provided  for  most  theorems,  since  they  can  be  found  in  any  text 
containing  Markov  processes,  such  as  Cinlar  [18]  or  Chung  [17]. 


Def.  2.1:  The  real-valued  stochastic  process  {Xt>  t  >_  0}  is  a 

continuous-time  Markov  Process  iff  Vn  and  ¥t  <  t?  <...  <  t 
-  i  z  n 

<  Vl' 


Be8(R),  the  linear  Borel  sets. 

This  Is  often  stated  as  the  future,  conditioned  on  the  present,  being 
Independent  of  the  past.  We  remark  that  in  the  above  definition  Xt 
may  be  a  vector,  which  is  in  fact  the  case  that  we  will  consider  later. 


Def.  2.2:  A  function  p  from  R+  xR  x8(R)  into  [0,1]  is  called  a 
transition  probability  function  (tfif)  iff 

(1)  P (t,x,-)  is  a  probability  measure  on  8 (R)  ¥t,x 

(2)  For  each  BeB(R),  p(.,.,B)  is  product  measurable  with 
respect  to  8(R+)  x  B(R) 

(3)  P  satisfies  the  Chapman-Kolmogorov  relation: 

VfieB(R),  s,t  R+, 

P(t+s,x,B)  »  /  P(s,x,dy)p(t,y,B) 
yeR 
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A  tpf  gives  the  probability  of  transitioning  from  some  starting 
position  to  some  Borel  set  in  some  fixed  time. 

Def.  2,3:  We  define  the  Markov  Process  to  have  stationary  transition 

probabi 1 ities  iff 


P[Xt+s  =BjXt]  =  p(s,  Xt,  B)  a.s. 

To  construct  a  Markov  Process  with  stationary  transition  pro¬ 
babilities  all  that  is  needed  is  the  tpf  and  an  initial  distribution 
measure  *(•)  on  SIR)  as  follows  from  the  following  theorem. 


Theorem  2.1:  Given  a  transition  function  p  ,  define  the  following 

finite-dimensional  distribution  functions  for  0  =  t  <  ...<  t  : 

o  n 


=  /  •  •  •  I  *(dx JpUi-t  ,x  ,dx1)p(t^-t1,x1,dx2) 

B  x . . . xBn  0  00 

o  n 


where  n  is  some  initial  probability  measure  on  6(R).  Then 

(Ft  }  is  a  consistent  family  and  hence  by  the  Kolmogorov 

Lo  n 

Consistency  Theorem  we  are  guaranteed  a  process  { Xt >  on 
n  =  {Xt>  has  stationary  transition  probabilities  and  is 

a  Markov  process. 

For  convenience,  we  will  label  P^  the  measure  constructed 
on  (r[0*°°) f  8(R^’”) )  from  *  and  p(t,x,B)  .  When  w(B )  =  <5X(6) 

then  we  write  P„  =  P  . 

x  ■* 
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2.2  Strong  Markov  Processes 

Suppose  the  Markov  process  is  defined  on  a  background  probability 
space  (£J,8,P).  Suppose  that  on  this  probability  space  we  have  an 
increasing  continuum  of  Borel  fields  8^.,  t  >_  0  ,  so  that  8tc8t 
if  s  >  t. 

Def.  2.4:  A  random  variable  T:i>+[0,«>]  is  a  stopping  time  with 
respect  to  {8t>  iff  (T  <  t)6  8 ■  The  pre~T 

sigma  field  is  defined  by:  8j={Ae8|a(T  <  tJeS^t}. 

Intuitively,  Bj  encompasses  all  of  the  information  obtained  by 
observing  the  process  up  until  time  T  .  It  can  be  easily 
checked  that  8y  is  indeed  a  Borel  field.  The  previous 
definitions  have  been  independent  of  the  Markov  process.  We 
now  suppose  that  £ 8^ }  is  the  set  of  Borel  fields  generated  by 
the  Markov  process,  i.e., 

8t  =  8(xs»  s  1  t),  where  8(XS,  s  <_  t)  is 

the  smallest  Borel  field  generated  by  the  random  variables 
{ X$ ,  s  <_  t}.  With  this  in  mind,  we  will  define  a  strong 
Markov  process: 

Def.  2.5:  A  Markov  process  {Xt>  is  said  to  be  strong  Marxov  if  V- 
stopping  times  T  and  ¥- xeR,  V- AeS  (R^0’°°^),  P_[X(*  +T)eA| 8j]  * 
Px(t)[x(-)£A]  as  Px  . 
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Thus  we  see  that  a  strong  Markov  process  is  one  for  which  the 
Markov  property  holds  for  stopping  times. 

2.3  Finite  State  Markov  Chains 

When  the  state  space  is  a  subset  J  of  the  integers,  we  call  the 
process  a  Markov  Chain  and  we  can  package  the  tpf  into  a  convenient 
matrix  form,  known  as  the  transition  matrix,  defined  as  follows: 

Def.  2.6:  A  function  P(t)  =  (P^-ft))^  .£j  is  called  a  transition 
matrix  if 

(i)  P^-U)  >0,1  P1j(t)  <  1  V-t  >  0 
J 

(i i )  P . . (t)  is  measurable 

^  j 

(iii)  The  Chapman-Kolmogorov  relation  holds: 

P(t+s)  =  P(t)  P(s) 

If  in  (i)  we  have  I  Pii(t)  <  1  we  ca^  p  substochastic.  If, 
j  J  ' 

moreover, 

(iv)  P(t)  •+•  I  as  t  \  0  we  call  P  standard. 

The  subsequent  discussion  will  only  involve  finite  state, 
standard,  stochastic  transition  matrices. 

We  have  seen  that  given  an  initial  distribution  and  a  transition 
matrix,  the  Markov  chain  is  probabilistically  determined  in  the  sense 
that  all  of  the  finite  dimensional  distributions  are  determined.  Since 
this  is  a  determining  class  (Billingsley  [9]),  then  two  different 
Markov  chains  with  the  same  initial  distributions  and  transition 
matrices  are  probabil istically  indistinguishable. 

From  the  statement  of  our  storage  problem,  however,  it  will  be  more 


natural  to  specify  the  infinitesimal  transition  probabil ities,  that  is, 
P.j j ( 6t )  for  fit  =  0  .  Specification  of  these  infinitesimal 
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probabilities  lead,  in  the  finite  state  space  case,  to  a  direct 
solution  for  the  transition  matrix  which,  along  with  the  initial 
distribution,  will  specify  the  Markov  chain. 

Def.  2.7:  The  waiting  time  random  variable  Wt  is  defined  by 

Wt  =  inf  (s  >  0|  Xt+S  t  Xt)  ,  and  Wt  =  «  if  the  set  is  empty. 
W  is  the  time  spent  in  state  before  jumping  to  a 
different  state. 


Def.  2.8:  The  jump  times  {T.  are  defined  by 
T0  =  0 

n  ■  «„ 

V,  •  X,  *  «T„ 

Def.  2.9:  The  sequence  of  states  visited,  {Xn}n«Q  are  defined  as 

*„  '  *(T„> 

{ Xn }  is  also  known  as  the  embedded  jump  chain. 


We  now  state  some  well-known  results. 


Theorem  2.2:  The  Markov  property  implies  that  Y-ieJ  3  x -e[0,»)  * 

-A.u 

P[Wt  >u|X^.ai]  =  e  ,u>^0 

Hence  the  holding  times  are  exponential  with  a  parameter  dependent 
on  the  state.  The  following  holds  for  joint  distributions: 


Theorem  2.3:  Vn ,  j  e  J ,  u  >_  0 

P  £Xn+i  =  J’»  Tn+i  '  Tn  ?  u  I  xo» •  •  •  »xn*To. .  •  •  .TR] 


-v 


Xn-J 


where  w  is  a  stochastic  matrix  3  i  0  ,  3  0  ,  it  1  *  1  .  By 

taking  conditional  expectations  of  both  sides  of  the  equation  with 
respect  to  X0,  ....  Xn  we  can  show  the  following: 


n  "xx.  ,uj 

Theorem  2.4:  P[T.  -  T.  >  u,,  j=l,...,  n|X0,...,  Xn]  =  n  e  j-1 
j  J  j=l 

Hence  the  times  between  jumps  (T^.-T^  }j<n  are  condit,ional  ^ 
independent  and  exponential 1y  distributed. 

As  a  corollary  of  this  theorem,  note  that  if  the  sequence  of  states 
visited  is  deterministic,  then  the  times  between  jumps  are  uncondition¬ 
ally  independent. 

By  conditioning  on  the  time  of  the  first  jump,  we  obtain  the 
Kolmogorov  Backward  Equation: 

Theorem  2.5:  V-i ,  j ,  t  _>  0 

Pii(t)  =  6ij  e  +  /xie  1  l  wikPkj(t's)ds 
1J  1J  o  keJ  1 

fu*  j 

By  conditioning  on  the  time  of  the  last  jump  before  t  ,  we  obtain 
the  Kolmogorov  Forward  Equation: 

Theorem  2.6:  ¥-1 ,  j ,  t  >_  0 

-At  t  “A . ( t-s ) 

ptj(t)  -  5ije  *  i  J,  pik(s>Vs*«« 

Q 

k/j 

Although  the  forward  equation  can  lead  to  difficulties  with 
explosive  processes,  that  is,  processes  for  which  P[XteJ]  <  1  for 
some  t  ,  finite  state  Markov  chains  are  non-explosive  and  if  the  jump 


matrix  is  irreducible  both  the  forward  and  backward  equation  possess 
the  same  unique  solution. 
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Theorem  2.7:  ¥-i,j,P_(t) 


has  a  continuous  derivative  and 


P'(0)  =  Q 


P'(t)  =  QP(t) 

l  -xi ,  i=j 

where  Q. . .  = 

U  |  Af*^. 


n  J*  j 


from  which, 


Def.  2.10:  The  matrix  Q  is  called  the  generator  of  the  Markov  chain. 

-1  t-h 

The  parameters  give  the  mean  holding  times  in  the  i  state  and 
n.j j  gives  the  probability  of  a  jump  from  state  i  to  state  j  in  the 
jump  chain. 

We  note  that  any  matrix  Q  such  that  Q  1  =  0  ,  q..  >0,  i/j,  and 

■**  **  I J 

<0  is  a  generator  for  a  Markov  chain. 

In  the  subsequent  analysis,  we  will  proceed  by  starting  with  a 
generator  for  the  Markov  chain  net  input  process  Xt  .  This  generator 
will  be  specified  by  the  infinitesimal  transition  probabilities.  We 
will  then  establish  results  for  the  contents  process  which  is  derived 
from  the  net  input  process.  We  will  maintain  as  much  generality  as 
possible  throughout.  The  most  general  setting  possible  is  that  on  an 
arbitrary  number  of  states  with  a  generator  as  described  above.  We 
will  refer  to  this  as  the  General  Case.  In  many  instances  this  is 
much  too  general  to  obtain  explicit  results,  and  so  we  restrict  our¬ 
selves  to  a  Markov  chain  on  two  states,  which  we  call  the  Two-5tate  Case. 
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In  this  case  we  take  as  a  generator 


As  we  discussed  earlier,  the  generator  by  itself  does  not  specify 
the  Markov  chain;  we  also  need  the  initial  distribution  tt  .  In  many 
cases  we  will  assume  that  it  is  the  stationary  distribution,  which  we 
discuss  next. 

2.3.1  Stationary  Distribution 

From  the  Markov  property  it  follows  that  if 

^•(t)  =  P[Xt=i]  and 

!  "  ("i}i=i 

then  ir" ( t )  *  n'  (O)P(t)  . 

Now,  suppose  that  there  is  a  time-independent  solution,  say  it  , 
to 

X'  *  X'  P{t) 

then  if  tt  ( 0 )  =  tt,  it  is  clear  that  ir(t)  t  >0  . 

In  this  case  the  Markov  chain  is  in  probabilistic  equilibrium  which 
we  call  (strict)  stationarity.  and  tt  is  known  as  the  stationary 
distribution. 

The  condition  for  existence  of  a  solution  of  the  system  above  can 
be  related  to  the  generator  by  the  following: 

Theorem  2.8:  x'P(t)  =  x'H  iff  x'Q  =  0  .  If  Xn  is  irredu¬ 

cible,  the  solution  is  unique. 

Hence  the  stationary  distribution  is  the  normalized  left  eigen¬ 
vector  corresponding  to  the  zero  eigenvalue,'  whose  right  eigenvector 
is  1  . 
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Now,  suppose  that  lim  P,,(t)  *  *.  <  ,  independent  of  i  .  Then, 

t-KX>  J 

in  matrix  notation, 

lim  P(t)  =  In' 

t-«o 

and  lim  n'(t)  =  rr'(O)  lim  P(t)  *  -"(0)  W'  = 

**  t  “♦<K> 

We  call  it  the  limiting  distribution  of  the  chain,  and  it  is  ree 
of  the  initial  distribution.  We  will  show  in  the  next  section  that  for 
our  case  the  limiting  distribution  is  identical  with  the  stationary 
distribution. 

Of  interest  is  the  rate  of  convergence  to  the  limiting  distribution 
for  the  net  input  process.  We  discuss  this  next. 

2.3.2  Rate  of  Convergence  to  the  Limiting  Distribution 

The  computation  of  limiting  distributions  and  results  on  the 
rate  of  convergence  are  most  easily  handled  by  algebraic  methods.  An 
excellent  reference  for  algebraic  methods  in  Markov  chains  is  Karlin 
[28].  The  basic  tool  used  is  the  spectral  decomposition  theorem  for 
matrices. 

Theorem  2.9:  If  Q  is  an  nxn  matrix  with  distinct  eigenvalues 

6 .  ,  i=l,...,n  and  right  eigenvectors  t^  ,  i=l,...,n  then  Q 
admits  the  spectral  decomposition 


r-l 


Q  =  TAT 

where  T  *  (t1,...,tn)  and  a*  diag  (e^). 

Now,  the  time-dependent  transition  matrix  was  given  by 


P(t)  =  e 


=  .Qt  .  V  fiV 


i=Q 


i! 


Substituting  Q  *  TAT"1  into  the  above,  we  obtain, assuming  distinct  s^s, 

•  i  ♦  t 

1*1  U 


I  +  T(eAt-I)T‘l 


T  AtT-l 

*  Te  T 
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9  .  t 


where  eAt  =  diag  \e  1  J . 
Hence 

P(t)  =  reAtT'1  . 

Now,  writing  T"1  = 


Tn 


we  can  write  the  above  as 
n 


p(t)  *  j,  *i  Z\ 


e,t 


and,  since  the  stationary  distribution,  tt  ,  is  the  left  eigenvector  of  Q 
corresponding  to  the  zero  eigenvalue  and  right  eigenvector  1  ,  we  have 

n  e.t 

P(t)  *  1  w'  +  I  t.  r;  e  1  . 

1-2  1  1 


We  now  prove  the  following  lemma  concerning  the  nonzero  eigenvalues  of 
a  generator. 

Lemma  2.1:  The  non-zero  eigenvalues  of  a  generator  have  negative  real 
part. 

Proof:  If  9  is  an  eigenvalue  then  0  satisfies  Qx  =  ex 
for  some  x  t  0  .  This  says  that  if  Q.  =  \  q.., 

1  m  1J 

n 

qi'JXJ  =  9xi  •  n 


or  -Qixi  +  l  q..x  =  9x . 
11  j^i  1J  J  1 


so  that  l  q.  x .  =  (e+Q.)x. 

j^i  J  11 

Now,  let  xk  *  V | x j  | 

l<j<n 
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Then,  for  i=k  we  get 


rx . 


|X, 


(e+Qk)  *  l  q,s.  (71- ]  with  \-M  <_  1 

j^k  \xk )  Xk 

|X,, 


so  then  |e  +  QJ  <  l  \%A  IrH  <.  >  ;qki| 

j>k  KJ  k  jh 


■  C, 


Hence  |8+  Qk|  1  Qk  »  or  |S  ■  (~Qk) l  1  Qk 
This  says  geometrically  that  8  lies  within  a  circle 
centered  at  -Qk  with  radius  Qk  on  the  complex 
plane,  so  that  Re(e)  <  0  . 


With  the  help  of  the  lemma  we  see  that 

Pit)  ■  H-  *  f  t(r;  e'1'1  ^ 

i=2  '  1  - 


with  the  rate  of  convergence  being  exponential  and  governed  by  the 
largest  non-zero  eigenvalue. 

Now,  if  ir'(o)  =  n'  ,  then 

Qt  r  oV 

n'(t)  *  ir^P(t)  =  ir'e  *  if'  l 

-  1-0  • 


I  +£  jt'QQ 
i=l  ' 


i-i  t‘ 


*  ir'I  +  0  * 


since  »'  Q  *  0 

So  we  see  that  the  limiting  distribution  w  satisfies  x'P(t)  =  x' 
and,  by  uniqueness  of  solutions  in  the  irreducible  case,  coincides 
with  the  stationary  distribution. 
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2.4  Weak  Convergence 

We  present  here  a  brief  discussion  of  the  theory  of  weak  convergence 
which  will  be  used  subsequently.  For  a  thorough  treatment  of  the  subject, 
the  book  by  Billingsley  [9]  is  recommended. 

We  begin  with  a  complete  separable  metric  space  (S,p).  and  we 
suppose  that  S  is  the  Borel  c-field  generated  by  the  open  sets  under 
p.  The  concept  of  weak  convergence  involves  convergence  of  sequences 
of  measures  defined  on  the  metric  space  as  follows: 

co 

Def.  2.11:  Suppose  IP  }  are  probability  measures  defined  on  (S,p). 

n  n=0 

We  say  that  { Pn }  converges  weakly  to  PQ  ,  written  Pn=>  P  iff 

/ «  p„  ^  /  fd  po 

for  all  bounded  continuous  real  valued  functions  on  S  . 

Weak  convergence  of  stochastic  processes  is  defined  by  considering 
the  weak  convergence  of  the  induced  measures.  To  do  this,  it  is  of 
benefit  to  consider  random  elements. 

Def ._  2.12:  X  is  a  random  element  of  (S,S)  if  3  a  probability  space 
(fl.B.P)  such  that  X  is  a  measurable  map  from  (q,8,P)  into 
(S,S). 

With  this  definition  in  mind,  we  can  define  weak  convergence  of 
random  elements  by 

Def.  2 ■ 1 A :  If  {Xnl  is  a  sequence  of  random  elements  on  (S,S),  then 

X  converges  weakly  to  X.  iff  P*X':  -^•X"1,  where  P-X*1 
M  u  n  o  n 

represents  the  probabil ity  measure  on  (S,S)  induced  by  Xn  . 

This  type  of  framework  serves  for  any  dimensionality.  For  example, 
if  (S,p)  *  (R,  | x-y | ),  then  a  random  element  is  a  random  variable  in 


k  Z  ic 

one  dimension.  If  (S,p)  =  (R  ,  [  7  (x.-y.)  P),  then  a  random  element 
is  a  random  vector.  For  our  analysis  of  weak  convergence  of  the  contents 
process,  we  will  work  in  the  space  S  =  D[0,»)  ,  the  space  of  right- 
continuous  functions  on  [0,»)  with  finite  left  limits.  Before 
discussing  the  D[C ,<*»)  metric,  we  will  discuss  the  simpler  space 
D[0,1]  .  Under  the  uniform  convergence  metric. 


p*(x,y)  =  sup  I  x ( t ) -y ( t ) | 

0<t<l 

the  space  D[0,1]  is  not  separable.  A  metric  under  which  0[0,1]  is 
separable  is  the  Skorohod  metric, 


p(x,y)  =  inf  p*(x,e)Vp*(x,yx) 

XeA 

where  A  =  {x:[0,l]  ■*.  [0,1]  |  x(0)  =  0  ,  X(l)  =  1,  x  is  continuous, 
one  to  one,  onto,  and  strictly  increasing  }  and  e  is  the  identity 
map.  Hence  A  comprises  the  set  of  time  transformation  of  [0,1]. 

The  idea  behind  the  Skorohod  metric  is  to  make  functions  which  are 
'close'  after  a  sufficiently  small  time  transformation  also  close  in 
the  metric.  Unfortunately,  D[0,1]  is  not  complete  under  the  Skorohod 
metric.  However,  Billingsley  has  modified  the  Skorohod  metric  into  an 
equivalent  metric  under  which  D[0,1]  is  complete.  For  details  see 
Billingsley,  pg.  112-113.  Clearly,  the  same  development  holds  on 
D[0,k]  for  any  k  . 

For  processes  on  D[0,<»)  ,  we  desire  to  define  a  metric  so  that 
random  elements  XR  will  converge  weakly,  Xn  ->XQ  in  D[0,«)  iff 
Xn  “>x0  °£0.k]  for  any  k  • 

To  accomplish  this,  let  Ak  be  the  set  of  homeomorphisms  from 
[0,k]  onto  [0,k]  with  the  properties  described  earlier  for  A  . 
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For  x,ye0[0,k],  define  the  Skorohod  metric 


dk(x,y) 


inf  p* 
,  .  k 
XeAk 


(x,e)  Vp*  (x,y.x) 
k 


where  p* 

k 

Let  d 
(D[0,k],dk) 
so  that  if 
define 


is  the  uniform  metric  on  [0,k]. 
k  be  the  Billingsley  modification  to 
a  complete  separable  metric  space. 
xeD°°  then  x  =  (xlf  x2,  ...)  where 


dk  which  makes 

co 

Now  let  0  =  "  D.  , 
k=!  k 

xkeDk  .  For  x.ycO” 


V"'  -k  /  dk^xk,yk^ 

d  (x,y)  =  L  2  - — 


k=l 


1  *  VvV 


Then  (D"d  )  is  a  complete  separable  metric  space.  Now  define  the 

co 

projection  maps 


rk:D[0,  «)  -*■  Dk 

by  rR(x(t))=  x(t),  0  _<  t  _<  k  ,  and  let  be  defined  by 

'k(x)  =  (rk{x),  k  >_  1 }  . 

Then  ^(D)  is  a  closed  subspace  of  D°°  and  (t|«(D),  d»)  is  a 
complete  separable  metric  space.  To  finish,  define  d  on  D[0,»)  by 
d(x,y)  =  d^(x),ii>(y)). 

This  is  the  metric  that  makes  D[0,«)  a  complete  separable 
metric  space  and  under  which  convergence  is  equivalent  to  convergence 
on  0[0,k]  for  all  k  .  For  a  detailed  treatment  of  convergence  in 
this  space,  the  reader  is  referred  to  Lindvall  [29], 
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We  now  define  determining  and  convergence-determining  classes: 

Oef.  2.14:  UcS  is  a  determining  class  if  two  probability  measures  P 

and  Q  on  S  such  that  P  =  Q  on  u  implies  that  P  M  •  0 

is  a  convergence  determining  class  if  for  any  sequence  of  probability 

measures  on  S,{Pn}“ 
n  n=0 

P  [A]  -  P  [A]  for  all  AeUaP-LaA]  =  0 
n  o  o 

implies  that  P  — >P_  . 

r  n  o 

A  convergence  determining  class  is  always  a  determining  class.  In 
many  cases  the  reverse  is  also  true.  For  instance,  in  the  cases  of  R1 
and  Rk  discussed  above,  and  even  in  sequence  space  R°°,  the  finite 
dimensional  rectangles  are  both  determining  and  convergence  determining. 
Therefore  to  establish  weak  convergence  it  is  sufficient  to  establish 
convergence  of  the  finite  dimensional  distributions  which  can  be  done, 
for  example,  using  characteristic  functions.  The  essential  difficulty 
in  D[0,°o)  is  that  the  finite  dimensional  rectangles  are  not  convergence 
determining.  Therefore,  convergence  of  the  finite  dimensional  dis¬ 
tributions  is  necessary  but  not  sufficient  for  weak  convergence. 

Convergence  of  the  finite  dimensional  distributions  along  with 
the  notion  of  tightness  is  necessary  and  sufficient  for  weak  convergence 
in  D[0,“>)  and  most  of  the  classical  proofs  use  this  argument.  How¬ 
ever,  tightness  is  usually  quite  difficult  to  establish  in  specific 
cases.  In  our  proof  of  weak  convergence  for  the  process  in 
Chapter  VI  we  are  able  to  avoid  tightness  arguments  by  considering 
an  embedded  renewal  process  and  an  appeal  to  a  technique  similar  to 
that  used  by  Durrett  and  Resnick  [20]  in  discussing  weak  convergence  with 
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random  indices.  The  major  usefulness  of  weak  convergence  lies  in  the 
continuous  mapping  theorem,  which  we  now  describe. 

Theorem  2,10:  Suppose  that  X  ,  n  _>  0  are  random  elements  of  (S ,S) 
defined  on  (n,5,P).  Suppose  that  h:(S,S,p)  (S'.S'.p'),  i.e., 

h  is  a  measurable  map  from  S  into  another  metric  space  S'. 

Let  Disc  h  =  {xeS|h  is  discontinuous  at  x}  .  If  P[XQeDisc  h]  = 

0  and  Xn“>  Xq  than  h(Xn)—>  h(XQ). 

By  considering  useful  maps  onto  other  metric  spaces,  then,  weak 
convergence  of  derived  processes  can  be  easily  obtained  from  the  basic 
weak  convergence.  For  example,  consider  the  mapping  from  D[0,»)  into 
D2[0,°°)  defined  by 

ht(x)  =(  Vx(u),  ax(u)  ]  . 

V 0<u<t  0<u<t  / 

This  is  a  continuous  mapping  from  D[0,°>)  into  D2[0,°°),  so  that  if  we  can 
establish  that  xn(.)  ->xQ(*)  in  D[0,=°),  then  it  follows  that 

/vx0(u),  AxQ(u) 

\0<_u<t  0<u<t 

Now  consider  the  continuous  mapping  from  D2[0,“)  into  D[0,»>)  defined  by 
h(x,y)  =  x-y,  and  we  can  establish  that 

Vxn(u)  -  Axr(u)  VxQ(u)  -  AxQ(u) 
o<u<_t  0<u<_t  0<_u<_t  0<u<t 

so  that  the  range  function  of  the  random  elements  *n(*)  will  converge 
weakly  to  the  range  function  of  the  limiting  random  element  xQ(*)  . 


Vxn(u),  Axn(u)' 


0<u<t 


0<u<t , 


CHAPTER  III 


ANALYSES  OF  MOMENTS  AND  AUTOCORRELATION  FUNCTION 
FOR  THE  DOUBLY  INFINITE  DAM 

The  doubly  infinite  dam  is  of  interest  for  preliminary  sizing 
studies,  as  has  been  discussed  earlier.  Before  any  useful  information 
can  be  obtained  from  an  analysis  of,  say,  the  range,  a  particular 
model  must  be  entertained.  Not  only  must  one  model  be  selected  out 
of  various  competing  models,  but  also  once  a  model  is  selected  there 
may  be  several  parameters  which  must  be  estimated  in  some  fashion.  For 
example,  the  Markovian  models  considered  in  this  investigation  contain, 
in  the  general  formulation,  the  following  parameters: 
m  states  uiv**.,u^ 

(m-l)m  jump  probabilities  . 
m  mean  holding  times  1  . 

This  is  a  total  of  m(m+l)  parameters  which,  if  m  is  only 
slightly  large  makes  the  model  difficult  to  work  with.  The  two-state 
case  contains  two  rates  and  two  mean  holding  times  for  a  total  of  four 
parameters,  and  the  symmetric  case  reduces  to  two  parameters.  Although 
these  special  cases  undoubtedly  oversimplify  the  situation,  neverthe¬ 
less  they  contain  a  manageable  number  of  parameters  from  which  exact 
expressions  of  quantities  of  interest  may  be  explicitly  obtained. 

In  order  to  both  check  the  fit  of  the  model  and  perform  an  initial 
estimate  of  the  parameters  it  is  desirable  to  know  the  theoretical 
moments,  or  in  this  case,  moment  functions  for  the  model.  In  particular 
the  autocorrelation  function  is  of  benefit  in  the  fitting  of  stochastic 
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models,  since  it  can  be  matched  with  the  sample  autocorrelation  function 
of  the  data,  which  can  be  easily  measured. 

In  this  chapter  we  will  investigate  the  following  functions  whicr 
we  now  define,  for  the  general  case  (as  far  as  possible),  and  the  two- 
state  case. 

For  a  stochastic  process  Z^,  t  ^0  with  finite  second  moments, 
we  define  as  usual  the 

mean  function  mt  =  /uP[Ztedu] 

(2) 

second  raw  moment  function  m£  '  =  /u2P[Ztedu] 

/  >  \ 

variance  function  -  m2 

cross-product  function  mst;  =  //uoP[Ztedu ,Zsedv] 
covariance  function  Kt  s  s  mst  -  msmt 
autocorrelation  function  p*  _  s  K* 

-  .  t  >5  l )  j  l  S 

If  the  process  Zt  is  strictly  stationary  with  finite  variance, 

then  clearly  m,.  and  a2  will  be  constant  and  pf  will  be  a  function 

of  t-s  only.  If,  on  the  other  hand,  these  three  conditions  are  met 

then  Zt  is  called  second-order  stationary. 

We  will  now  investigate  as  much  as  possible  the  above  functions 

for  the  marginal  process  on  the  doubly  infinite  dam:  Xt,  the  Markov 
t 

chain,  and  C t  =  /  Xudu  ,  the  unrestricted  contents  process. 


3.1  General  Case 

3.1.1  Markov  Chain 


Suppose  that  the  Markov  chain  Xt  has,  as  generator, 

i ,  j=l , . . .  ,m 

\  'J 


Vij.  w 


and  suppose  that  it  is  defined  on  states  uj,...,um  . 


Let 
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ir.(t)  =  P[Xt=i]  as  before.  In  the  subsequent  formulae,  we  will  drop 
the  upper  and  lower  indices  on  the  summands,  which  will  always  be  m 
and  1  respectively.  We  assume  throughout  that  the  eigenvalues  of 
Q  are  distinct. 


mt  =  £  Vi  (t) 


"t.t+s  l  |  uiujP^Xt  =  1  ’  Xt<s  =  ^ 


(3.1) 


and 


=  l  l  u.u.P.^sJnJt)  by  the  Markov  property  (3.2) 
i  j  1  J  1J  1 


=  [  l  uJui-uJVtJir.U)  (3.3) 


For  the  covariance  function  we  obtain,  using  (3.1)  and  (3.2) 

Kt,t+s  =  l  l  MiwjPij(s),ri(t)  '  l  Vi(tJ  I  wj7rj(t+s)  (3>4) 
1  J  ^  J 

giving  for  the  autocorrelation  function 

l  I  UiUjP1j(s)irf  (t)  -  ^  u^U)  l  lijiTjU+s) _ 

U-j  “Uj  (  t+S  )lTj  (t+S  )J 
■  ,  J 


t,t+s 


(3.5) 


Now,  when  the  Markov  chain  has  its  stationary  distribution  r  ,  we 
see  that  ot  t+s  is  a  function  of  s  alone,  given  by 


l  l  VjPij(s)”i  -  (£Vi) 


(3.6) 


y  ui (ui'uj )iri7rj 
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Using  the  expression  for  P..(s)  from  section  2.3.2  this  reduces 


to 


V 


ps  ’ 


SEE  tkirkje 

i  t  j 


(3.7; 


or 


m  9  s 

p  =  2  Cke 

s  k=2  K 


(3.8) 


where 

_  l  1  tkirki'i“t"j 

k  y  “)<“r“j>vj 

Since  Re(ek)  <  0  ,  this  shows  that  p$-»0  as  s-*».  This  also 
shows  that  the  rate  of  convergence  is  governed  by  the  eigenvalue  with 
largest  real  part. 

3.1.2  Contents  Process 

t 

We  now  consider  the  contents  process  =  j  X^du  .  For 
this  analysis  we  also  assume  that  X^.  is  stationary.  is  the  net 
input  (input-release)  in  (0,t)  for  a  doubly-infinite  dam  whose  net 
rate  of  change  of  level  at  time  u  is  Xu  .  This  is  the  continuous 
time  analog  of  the  cumulative  sums  of  a  Markov  chain  (Xn,  n=0,l,2,...} 
as  considered  by  Odoom  and  Lloyd  [36],  Ali  Khan  and  Gani  [1]  and  others. 
Since  the  typical  procedure  is  to  observe  the  input  at  a  discrete  set 
of  times  {0,h,2h,...}  it  is  useful  to  know,  for  modelling  purposes, 
the  mean  and  covariance  function  of  the  increment  process 
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AhU)  = p  411  X(u)  du  =  C  (c  +  |)  -  C  (t  -  |)  (3.9) 

t-%h 

These  are  the  same  increments  treated  by  McNeil  for  a  special 
two-state  case  where  one  of  the  states  is  zero.  We  will  examine  A^(t) 
more  generally.  The  relevance  of  a  (t)  is  that  by  looking  at 
°t  t+h  ^or  t*ie  Ah^  variable,  we  will  be  examining  the  correlation 
between  adjacent  increments.  This  correlation  will  then  provide  an 
indication  of  the  degree  of  approximation  to  our  process  by  one  with 
independent  increments  for  which,  of  course,  t+h  =  0  ‘  For  the 
rest  of  this  section,  the  moment  functions  which  we  develop  will  be  for 
the  increments  of  the  process,  which  we  will  denote  by  an  argument  of 
h  in  the  function. 

For  the  mean  function  we  have 

t+  h 

f t+  7 

EAh(t)  =  El  x (u )  du  =  l  u^h 


For  the  second  raw  moment  function,  we  have 


mt2)(h)  =  E  Ah(t) 


2 


x(u)  du 


E[x(u)xU)]du  d£ 


h 

2  E[x(u)x(5)]  du  d5 


(3.10) 
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I 


where  we  are  now  integrating  over  the  half-rectangle  in  which 


u  >  5  . 


Using  (3.2)  and  the  results  of  section  2.3.2,  we  can  write  (3.10) 


h  ^*^h 


vj 


t_  2  ? 


ek(u-e) 


,hiru*  Jwidu  dc 

(3.11) 


After  performing  the  indicated  integrations,  (3.11)  becomes 


42)oo  »  2  £  £  Z 

^  _•  a  i _ o 


^  Vr„i  r  9kh 

>  ki  kj 


i  j  k=2  8 


+  l  l  lVjiriTrj  h2 

i  j 

For  the  variance  function  we  obtain 


e  ekh  -  1 


(3.12) 


a|(h)  =  m^  ^ 


Ex(u)  du 


■  *(21  -  [i  «,*i» ) 2 


(3.13) 


and,  substituting  from  (3.12),  we  obtain 


at(h)  -  2 


m  t  r  r 6 

i  j  k=2  0,,  L 


ekh  -  lj 


(3.14) 


i 


! 
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For  the  covariance  function  we  have 

r 


Kt,t+h(h)  *  Cov 


f  ?  rt+  j- 
/  x(u)  du,/  x(u)du 
.  ;  h  .  h 

L  ?  u  +  2 


t+  3h 

*/V  / 

L  +  7T 


‘♦r 


7  *-f 


Cov[x(u)x(?)]  dc  du 


and,  substituting  for  Cov[*,.]  from  the  numerator  of  (3.7), 


we  have 


Kft+h(h) =  y 

t  - 


t  +  £ 
c  ♦  2 


3h 


rltr  "  efc(s-“) 

/  h  2^  2-A^  e 
•/  +  h  i  j  k=2  KJ 


t  -  ? 


VjV" 


(3.15) 


The  integration  in  (3.15)  is  straight-forward  and  we  obtain  finally, 

h) 


-e 


Expanding  the  product  makes  the  dependence  on  t  disappear,  and 
we  have 


1  J  K=2  9k 


2ekh  e,  h 
1  +  e  k  -  2e 


J  uiujwi.  (3.16) 


Since  the  variance  of  the  increments  is  free  of  t  ,  we  can  use 
equations  (3.14)  and  (3.16)  to  write  the  correlation  function  as 
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McNeil  notes  that  for  his  two-state  model  with  ui  =  a  and  u2  =  °» 
pt  t+h  is  1'ree  a  '  ^his  becomes  obvious  from  (3.17),  since 
in  this  case  there  is  only  one  summand  and  the  term  becomes  a2  , 

which  cancels  from  the  numerator  and  denominator. 

We  close  this  section  by  checking  (3.17)  for  the  McNeil  model.  For 
the  McNeil  model  we  have  m=2,  ui=a  and  p2=o  ,  whereas  for  a  two- 
state  model  we  have,  as  will  be  derived  in  the  following  section, 
ttj  =  p/A+p  and  e2  =  -(x+p).  Hence 


the  same  as  found  by  McNeil. 

3.2  Two-State  Case 

In  the  two-state  case,  we  can  obtain  some  simple  expressions  for 
the  moment  functions. 

3.2.1  Markov  Chain 

In  this  case  the  generator  matrix  is 


which  admits  the  spectral  decomposition 
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=  (x+p)  ^owJ+Xuj) 


3  pX(X+p)  2 ( ui -U2 V 


mt,t+s  =  (x+o)‘ 


<s  =  px(x+p)  { u i _U2 )  e 


(pU  1+Xy2  )  +  Xp  (u  1"U2  )  e 

2  -(X+p)s 


•(x+p)s 


-(x+p )s 


Ps  -  e 


We  note  that  all  limiting  values  have  an  exponential  rate  of  conver¬ 
gence  with  (x+p)  as  a  multiplier  of  t  .  Thus  X  and  p  both 
contribute  equally  to  the  order  of  convergence.  In  addition,  note 
that  ps  is  free  of  u\  and  u2  . 

3.2.2  Contents  Process 

For  the  contents  process,  when  Xt  is  stationary,  we 

obtain 


m 


t 


(p+X)  1 


(py 1+Xy2 )  t 


(2)  _  r  r  ft 

t  “  l  l  ^iVi  l 


J5  Pii(s)ds  +  Pii(s)ds 


dC 


The  integrations  are  straight-forward,  and  we  obtain  finally 

(?)  7  -4  z\  -(x+p)t  1 

y  ;  3  -  2 (x+p )  px(m-y2)  [l-e  -(x+p)tj 


from  which 


o|  =  2 (X+p )  4pX(u1-y2)i 


I*  ”(x+p  )t  1 

Le  +  (x+p )t  -  1J 
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We  notice  from  the  above  that  the  variance  function  is  proportional  to 
the  square  of  the  difference  between  the  two  states,  and  the  large-t 
growth  is  approximately  linear. 

Rather  than  examining  the  correlation  between  arbitrary  points  of 
time,  we  again  consider  the  correlation  between  increments, which  is 
easily  obtainable  from  eq.  (3.17) 


t,t+h 


(h)  . 


2e->h  e2hl  t,-r0. 

1  ♦«  ■  -2«  l 

L _ J  i  .1  1  j  1 


S-h  t?r 

2  -02h+e  2  -1  ll  Ar1 

J  i  j  e2 


uiuj7Ti 


-(x+p)h 


2  (x+p )h  -  1  +  e 


-  X+p  hi  • 


the  same  as  found  by  McNeil.  Note  that  pt  t+h  is  free  of  the  states. 
This  is  a  slightly  more  general  result  than  the  one  reported  by  McNeil 
where  only  one  state  was  arbitrary. 


CHAPTER  IV 


ANALYSIS  OF  THE  RANGE  IN  THE  DOUBLY  INFINITE  DAM 

In  this  chapter  we  begin  the  range  analysis  for  the  unrestricted 
t 

contents  process  Ct  =  /  X  du  in  the  doubly  infinite  dam.  Let  T 

o  . 

be  an  exponentially  distributed  random  variable  with  mean  s-1  an 

independent  of  the  process  {Xt }  ,  and  let 


M$  =  supIC^ ,  0  £  u  £  T}, 


(4.1) 


m$  =  inf{Cu ,  0  £  u  £  T} 


where  as  before  { X t >  is  a  finite  Markov  chain  with  generator  Q  and 

Ct  •  /*  xudU  . 

We  shall  determine  the  joint  distribution  of  Mg  and  which 
in  principle  determines  the  joint  distribution  of 


Mt  =  sup{Cu>  0  £  u  £  t}. 


=  inftC^,  0  £  u  £  t>. 


and  hence  of  the  range 


A  closed  form  expression  for  the  distribution  of  R^  is 
difficult  to  obtain;  however,  we  shall  find  an  explicit  expression  for 
ER^  in  the  symmetric  case  and  investigate  its  asymptotic  behavior  for 
large  t. 


43 


44 


4.1  Joint  Distribution  of  H  ,  in 
_ s  s 

Define 

s(x,ui)  =  P[x+M$  <_  a,  x+ras  >_  0 1 XQ  *  u^]  ,  0  <  x  <  a 

(4.5) 


Then,  appealing  to  the  strong  Markov  property  of  (Xt,Ct)  as  shown  by 
Erickson  [20-1]  we  can  write  the  first  jump  equation  for  0  as 


e(x,u1-)  =  0  ( x ,  u  j ^  ^"e(x+}i.jU,u.j  )1T1-je 

J  *  .  . .  -  \ 


V 


x  .du 

i 


(0<x+u.u<a ) 


(4.6) 


where  p(x,u^)  represents  the  case  in  which  there  are  no  jumps.  For 
p(x,Ui)  ,  therefore,  we  can  write 


Hence 


pfx.Uj ) 


(^  <  0) 


(ui  >  0) 


p(x,Ui  ) 


Vs 


x.+s 

1 


1 

1 


e 


e 


( x 1 +s )  x/p. 


(u-j  <  0) 
(ui  >  0) 


(4.7) 


(4.8) 


(4.9) 

(4.10) 


9  is  a  bounded,  measurable  function  on  the  state  space  of  the 
bivariate  process  (Xt,Ct).  Hence,  as  shown  by  Brockwell  [12]  the 
system  of  integral  equations  possess  a  unique  bounded  solution  and 
9(x,p.j)  satisfies  the  Kolmogorov  backward  equations. 
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Writi ng, 

^(x.s)  =  1  -  9(x,ui)  =  1  -  P[MS  <_  a-x,  ms  >_  -x|xQ  =  ui] 

(4.11) 


we  find  by  differentiating  (4.6)  that 

—-^(x.s)  =  D" 1  (Q-s I )  i{/(x , s )  (0  <  x  <  a) 

where  «r(x,s)  =  ( (x,s) , . . . ,  i|/m(x,s)j 
and  D  =  d i a g ( u ) 
with  boundary  conditions 

H/i(a,s)  =  1  (ui  >  0) 

^■(O.s)  =  1  (u^  <  0)  . 

The  system  (4.12)  has  the  solution 

tl»(x,s)  =  exp[-D"2  (Q-sI)xj  i|<(0,s) 


(4.12) 


(4.13) 


4.2  Explicit  Solution  for  the  Two-State  Case 
Suppose  that  Xt  has  generator 


and  is  defined  on  the  states  m  >  0  and  u2  <  0  •  1°  order  to 

determine  exp[-0_1 (Q-sI )x]  ,  we  spectrally  decompose  M  =  -D-1(Q-sI). 
The  eigenvalues  are  given  by 

9  *  1  +  +  /  \  '(  e+s  +  x+s  V  s(A+p)  +  S2 

2  L  U2  WlJ  ~  /  4  y  u2  /  UlU2 
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We  note  that,  if  we  set 


then 

9=d+([d|  +  A ) ,  where  A  >  0  . 

Case  1 :  d  >  0  .  Then  9!  =  d  +  d  +  A  =  2d  +  A>0 

92=d-d-A  =  -A<0 

Case  2:  d  <  0  .  Then  91=d-d+A=A>0 

e2=d-(-d)-A=2d-A<0 
So  that  in  either  case  we  have 


. 9 1  >  0  ,  02  <  0. 

The  right  and  left  eigenvectors  can  be  taken  to  be 


tf  -  (A.A+S-9^1 ) 
t£  =  ( A , A+s-92u i ) 

rl  =  (iJjP  >W2(^  +  S-01Ul  )) 
r2  =  (  Uj  P,M2(A+S-92U!  )) 

Mx 

We  now  note  that  we  can  write  e  v  ,  v  an  arbitrary  column  vector,  as: 

Mx  9ix  P2X 

e  v  =  t^r^v  e  +  t2r2v  e 


Hence 


where 


Al  _  tn 

=  Ff2’  i'e"’  ti2Ai  =  tnA2 


(4.14) 
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and 


Bl  ^21 

IT  =  ’  i,e-’  t22Bl  =  t2lBz  • 

With  this  in  mind,  we  can  write  the  system  (4.12)  as: 

OjX  9 2X 

Mx’s)  =  A:e  +  Bje 

BjX  e2x 

4»2  (x ,  s )  =  A2e  +  B2e 

with  auxiliary  conditions  for  Aj ,  A2,  Bj,  B2  given  by 

e^a  62a  \ 


Ate  +  Bje 

A2  +  B2 


=  1 
=  1 


from  the  boundary  conditions 


and 


( X+S-6 i u i )A\  -  XA2 


from  (4.14) 


(x+s-02uj )Bi  =  XB2  ] 

Solving  for  b.  and  B  •  we  obtain 
e2a 

xe  -  (x+s-e2ui) 

1  e^a  eTa 

e  (x-*-s-9lUl )  -  e  (x+s-e2Ul) 


(4.15) 


(4.16) 


9 13 

xe  -  ( A+S-9 1 Mi  ) 

^  9  j  a  9  2a 

e  ( X+S-02U i )  -  e  (x+s-9iui) 


X  +  S-9  2p i 

X 


Bl 


The  exact  joint  distribution  of  M$  and  mg  are  easily  obtained  from 
(4.15)  and  (4. If),  from  which  the  exact  distribution  of  Rs  =  Ms  *  ms  * 
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the  range  to  time  T  can  be  obtained  in  principle.  An  exact  solution, 
however,  is  analytically  intractable.  Hence,  rather  than  looking  for 
the  exact  distribution  of  R$  ,  we  will  look  instead  for  its  expecta¬ 
tion,  since  by  linearity  of  the  expectation  operator,  we  can  obtain 

this  using  knowledge  only  of  the  marginal  distribution  of  M  and  m  . 

s  s 

4.2.1  Marginal  Distribution  of 

In  the  definition  of  ,  letting  y  *  a  -  x  ,  we  have 

«t»i (y.s)  =  1  -  P[M$  iy,  ms  >y  -  a|XQ  =  u.]. 

Therefore,  the  marginal  distribution  of  M$  is  given  by 

P[Ms  >  y]  =  lim  ^(y.s)  . 

3-k» 

Since  ej  >  0  and  e2  <  0  ,  and  noting  from  4.16  that  lim  B.  exists 

a-*»  1 

and  is  finite,  i=l,2,  then  from  4. 15  it  follows  that 

6i(a-y) 

lim  ij,. (y.s)  =  lim  A^e  ,  i  =  1,2 

3-h»  d-KJo 

and,  after  evaluating  this  limit  we  have 

-9iy 

My. *)  =  e 

A+s-ejui  -Qjy 

My»s)  3  — i —  e 

Therefore,  if  the  initial  rate  is  positive  then  M$  is  exponentially 
distributed  with  parameter  Bi  .  If  the  initial  rate  is  negative,  then 
Ms  has  a  truncated  exponential  distribution  with  parameter  ei  and  a 
mass  at  0  of  size 
eiin-s 
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Therefore  it  follows  that 


E[Ms|u,]  =  ej1 

E[Ms|y2]  =  (01A)”1(A+S-0iu1 ) 

4.2.2  Expected  Value  of  R$ 

By  the  law  of  total  probability 

ERS  =  E(Ms-ms)  =  [e(Ms |UI )  -  E(m$ |ui )]  [p  X(o)  =  ux  ] 

+  [e(M$ fu2 )  -  E(«ns [u2 )]  [p  X(o)  =  U2  J 

However,  from  consideration  of  synmetry  we  have 

ECms i u 1 3  -  -  e£m$  | -W  X ] 
eOs|u2]  =  -  E[Hs|-u2] 

Hence 

E[Ms-ms]  =  [E(Ms  |ui )  +  E(Ms|-m)]  P[X(o)  =  Ul]  (4.17) 

+  [eCMs|u2)  +  E(MS | -U2 )]  P[X(o)  =  u2] 


Stationary  Case 

We  now  evaluate  ERS  when  X(t)  has  the  stationary  distribution, 
which  in  the  two-state  case  is  n'  =  ^p(a+p)  1  A(x+p)_1j. 

Substituting  the  appropriate  values  into  (4. ’7)  we  obtain: 


ER. 


(A+p) 


-1 


-1 

01 


(p+A+S-0 1 y i )  +  8i  ^A+pA  (A+S+01y2)j  (4.13) 


where 
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9i  * 


e 


1 


1  £±i  +  hll  +  / 1  |P±§.  +  AilV  .  lLit£  )+s2 

2  L  U2  Ui  J  /  4  l  U2  Ui  /  UJU2 


1  fill  +  x+s  +  7 1  ( £±£  +  A+i )  _  s(x+p )+s2 

2  L  PI  U2  ^  4  \  lii  U 2/  UlU2 


4.2.3  The  Laplace  Transfonn  of  ERt 

If  mR(t)  =  ERt  .  then  since 

»  -st 

ER$  =  /  mR(t)se  dt  , 

it  follows  that  the  Laplace  transfonn  with  respect  to  time  of  ERt.  say 

L 

A 

mR(s)  ,  is  given  by 

»R(s)  3  . 


Although  the  general  two-state  solution  is  easily  obtained  from  (4.18) 
we  will  further  restrict  our  attention  to  the  symmetric  case  to  perform 
an  exact  inversion: 

Symmetric  Case  ui=-U2=b,A=p=a 


In  this  case  we  obtain  after  simple  algebra  that  e 


l  =  =/- 


2ab+s^ 


and 


mR(s)  *  ba’1 


2a  +  s 
/2as+s2 


b2 
(4.19) 


We  can  write  mR(s)  as 

mR(s)  =  ba“ Xs~ 1  [(2as ) -Js(2a-»-s )  ( l+s/2a ) _is-l]  . 

Expanding  (l+s/2a)”Js  in  a  Taylor  series  about  o  ,  we  obtain 
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'  /  a  4  •  n ! 


Hence,  combining  like  powers,  we  have 


«•♦*>(' ♦&)J‘ ■*..(' 


+...+• 


(.1)0-1  l‘3...(2n-3) 

v  ;  Tn-nr^r  rr~ 

4  a  (n-1 ) ! 


+  (_1)n  2«a-l «3. . . (2n-1 ) 


4nannl 


s  +  ... 


Noting  that 


l/.lyM2n~3.? _ .  2*a,l  *3  ...  (2n-l )  2-l»3  ...  (2n-3) 

4n-lan-l(n_i)!  4nanni  ~  >tO»n-l., 


4  m  M  I 

a  n! 


4  a  '  n! 


we  have 


(2a+s)(l+s/2a)-*5  =  2a  ♦  1  s  +  Y  M)"'1  Ml3 ...(2n-3)  n 

/  .  '  '  n  r»  _  7  a 


n=2 


4nan“1  n! 


from  which  we  get 


1  -3/2  -3/2 


mR(s)  =  —  b  a  s  (2a+s)(l+s/2a)  -  ba‘ls"1 

VC 

=  /?  b  s’3/2  -  b  a"ls"1  t|ba  "3/2  s'*5 

oo 

*;STr:2  Z  1 3-(2j*i)  (i\^~ 

j=0  4J(j+2)!  \a/ 


(4.20) 


valid  for  0  <  |s|  <  2a  . 

4-2*4  Explicit  Inversion  of  Laplace  Transform 

Equation  (4.2.0)  gives  an  expansion  for  mR(s)  in  increasing 
powers  of  s.  Although  such  a  form  can  not  be  used  for  term  by  term 
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inversion,  it  will  be  useful  to  derive  an  asymptotic  expansion  which 
will  be  discussed  later.  We  now  perform  an  expansion  in  terms  of 
decreasing  powers  of  s  from  which  explicit  inversion  will  be  possible. 


Returning  to  (4.19)  we  can  write  it  as 


valid  for  |^-|  <  1  ,  i.e.,  |s|  >  2a. 

Equation  (4.21)  provides  a  power  series  expansion  for  m^(s)  which 
is  absolutely  convergent  in  a  neighborhood  of  infinity.  In  order  to 
perform  term-by-term  inversion,  we  appeal  to  a  theorem  found  in 
Ooetsch  [19]  which  we  now  state: 

Th:  When  f(s)  can  be  expanded  in  an  absolutely  convergent  series 
for  | s |  >  R  of  the  form 


where  the  x  form  an  arbitrary  increasing  sequence  of  numbers 
0  <  \Q  <  then  the  inversion  can  be  executed  term  by 

term: 
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The  series  for  f(t)  converges  for  all  real  and  complex 
t  t  0  . 


The  conditions  of  the  theorem  are  met  by  (4.21),  with 

an  =  (  n  )  2nan  and  *n  =  n  +  1  . 
Therefore,  we  have  for  the  inversion 

h  \ 


mR(a,b,t)  =  b  £  ( k  )  2a* 


k! 


(4.22) 


The  series  in  (4.22)  can  be  expressed  in  terms  of  a  confluent  hyper 
geometric  series  which  we  now  define: 

Def .  The  confluent  hypergeometric  series  with  parameters  cteR  and 
ygR  is  denoted  by  F(a,y,x)  and  is  given  by 


where  (a).  =  a(a+l ) . . . (a+j-1 )  for  j=l,2,... 

J 

=  1  for  j  =  0  . 

The  confluent  hypergeometric  series  converges  for  all  x  (Buchholz 
[15]),  and  has  been  tabulated  in  [27]. 

Returning  to  (4.22)  we  have 


«(a,b,t)  *  j  S  (k  ) 
t  a  k=l 


0k  k.k 
2  a  t 

k! 


Hence, 


t(a,b,t)  =  f  [F(-is,l  ,-2at)-l] 


(4.23) 


The  confluent  hypergeometric  function  for  a  negative  argument 
behaves  like  a  very  slowly  converging  alternating  series.  This  makes 
it  difficult  to  compute  numerically.  However,  we  can  use  Kummer's  first 
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formula  found  in  Buchholz  to  express  it  in  terms  of  a  hypergeometric 
function  with  positive  argument.  Kummer's  first  formula  is 

F(ct,e,x)  =  exF(e-a,S,-x)  . 

Using  this,  we  can  re-write  (4.23)  as 

b  r  _2at 

mR(a,b,t)  =  J  [a  F(3/2,l ,2at)-l  . 

For  numerical  calculation  on  a  computer  this  form  is  useful,  for 
the  confluent  hypergeometric  function  is  now  a  monotcnic  series  which 
converges  rather  quickly,  and  the  internal  computer  algorithms  can  be 
used  to  calculate  the  exponential  term. 

We  present  one  last  form  for  mR(a,b,t),  using  formula  (14),  page 
7,  of  Buchholz: 

-2at  *  2at 

mR(a,b;t)  =  — : -  /  ev  vi5(2at-v)'3/2  dv  (4.24) 

a^r(-Js)  { 

4.2.5  Asymptotic  Expansion  of  Laplace  Transform 

Equation  (4.23)  provides  the  exact  mean  of  the  range  to 
time  t  in  terms  of  a  known  function.  However, the  asymptotic  behavior 
is  not  evident  from  this  form.  To  study  the  asymptotic  behavior,  we 
make  use  of  the  technique  of  asymptotic  expansion  of  the  Laplace  Trans¬ 
form,  following  Doetsch. 

We  first  define  an  asymptotic  expansion: 

Def .  A  function  ${z)  is  said  to  have  an  asymptotic  expansion 

CO  00 

l  ij/  (z)  as  z-*»  ,  written  <f>(z)  =  l  ip  (z)  ,  iff 

v*0  v=0  v 

(j>{z)  -  l  I pjz)  *  o(<i,n(z))  as  z-«  for  every  n=l,2 . 

Another  way  of  stating  this  is  to  say  that 
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<p(z)  -  l  i j»v(z) 

- ^  -  0  as  z-  . 

n 

This  says  that  the  error  i>(z)  -  [  ip  (z)  not  only  goes  to  zero,  but  it 

v=o  v 

is  of  lower  order  of  magnitude  than  the  last  term  in  the  sum. 

We  now  state  a  theorem  found  in  Doetsch  which  can  be  used  to  obtain 

an  asymptotic  expansion  for  m^(t): 

Theorem:  If  f(s)  is  the  Laplace  transform  of  f(t)  ,  and  f"(s)  can 
be  expanded  in  a  neighborhood  of  aQ  in  an  absolutely  conver¬ 
gent  power  series  with  arbitrary  (perhaps  non-integral) 
exponents: 

oo  ^ 

f(s)  *  [  Cv(s-a0)  v  ,  -  N  <  <  ...-►« 

v=0 


(where  N  is  a  positive  integer)  then  the  following 
asymptotic  expansion  for  f(t)  is  valid  for  t-»>  : 


f(t)  =  e  0  22 

v=o  r(-xv) 

with  ~  }  =  0  if  xv  *  0,1,2... 


We  wish  to  apply  this  theorem  to  (4.20)  with  aQ  ■  0  ,  but  first  we 
must  remove  the  singularities  at  zero.  To  do  this,  define 


f(s) 


mR(s)  -  /2  b  a 


-*3  -3/2 


+  b 


-i-i  , 
a  s  -  -r-  b 
4 


a-3/2 


a 

f(s) 


/?ba'2 

16 


<30 

z 

ja0 


ll3,_.,12j^lj 

4j(j+2)! 


Then 
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In  particular,  since  t'^-^O,  then  as  t-*» 


Notice  that  if  b  =  vT  then 
mR^-/7th 


(4.27) 


which  is  the  value  for  the  Wtener  process.  This  is  far  from  coinciden 
tal:  In  Chapter  VI  we  will  prove  that  under  these  conditions  the 
process  will  in  fact  converge  weakly  to  the  Wiener  process. 

4.2.6  A  Note  on  the  Hurst  Phenomenon 

For  the  Wiener  process  it  is  well  known  that  the  maximum 
random  variable  has  distribution 

P[Mt  <_  x]  =  1  +  2  <f(0,t;x) ,  x  >_  0 

where  4(0, t;x)  is  the  normal  distribution  function  with  mean  0  and 
variance  t  . 

Hence  we  obtain  for  the  mean,  say  EM^ 

x2 

EH  •  f"  -2i-  e  2t  dx  =  /-  t1 

Since  the  expected  range  is  given  by 
ERB(t)  =  EMt  -  Emt 
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and  since  for  the  Wiener  process  Emt  -  -  EMt  ,  it  follows  that 

ERg(t)  =  2EMt 
or 

ESB(t)  '  A  tH 

By  appealing  to  the  continuous  mapping  theorem  one  can  hope 
that  this  well-known  result  for  the  Wiener  process  will  also  hold 
asymptotically  for  any  process  which  converges  weakly  to  standard 
Brownian  motion. 

H.  E.  Hurst  [26],  in  studying  large  amounts  of  data  involving 
streamflows,  rainfall,  pressures,  etc.,  found  that  ERn  appears  to 
vary  as  na  with  .69  <_  a  <_  .80.  The  average  values  of  a  observed 
by  Hurst  was  .72.  Since  his  discrete  time  models  for  these  processes 
were  assumed  to  be  made  up  of  sums  of  independent  radom  variable,  then 

5 

the  theoretical  behavior  of  ERn  should  have  been  as  n'  .  This 
discrepancy  became  known  as  the  Hurst  phenomenon  and  was  commented  on 
in  the  introduction. 

Although  no  completely  satisfactory  answer  has  been  provided  to  the 

phenomenon,  at  least  two  theories  have  been  suggested.  The  first  is 

that  the  assumption  of  independence  is  not  valid,  and  the  observed 

growth  of  ERn  is  the  result  of  serial  dependence  in  the  series  of 

summands.  However,  as  noted  by  Moran  [33],  such  a  dependence  would 

have  to  be  of  a  very  peculiar  kind,  since  with  all  plausible  models 

the  large  time  behavior  is  approximated  by  an  additive  process,  whose 
0  5 

growth  rate  is  n  . 

The  second  theory  is  that,  due  to  dependence  of  the  summands,  the 
series  studied  by  Hurst  are  not  of  sufficiently  long  duration  to  exhibit 


the  asymptotic  behavior,  and  hence  what  he  was  seeing  was  the  p re- 
asymptotic  behavior.  Several  models  have  been  suggested  which  exhibit 
this  behavior  in  discrete  time. 

In  our  continuous  time  model,  the  analog  of  dependence  of  the 
summands  is  dependence  of  the  increments  of  Ct  .  Asymptotical ly  we 
have  shown  that  ERt  behaves  like  j—  t  .  However,  it  is  of  interest 
to  Investigate  ERt  for  finite  t  to  see  if  the  Hurst  phenomenon 
appears  pre-asymptotically. 

In  Figure  4.1  is  a  plot  of  mR(t)  vs.  t  for  both  our  Markov 
model  with  a  =  b  =  1  and  for  Brownian  motion.  In  Figure  4.2  is  a 
plot  on  a  logarithmic  scale  of  mR(t)  ana  t*^  vs.  t  for 
•1  <_  t  <_  2  .  From  this  plot  it  appears  that  the  rate  of  growth  of 
mR(t)  is  closely  approximated  by  a  power  of  t  ,  and  within  this  range 
of  time,  the  power  of  t  is  close  to  Hurst's  average  of  0.72. 

The  rapid  convergence  to  the  asymptotic  value  is  apparently 
inconsistent  with  the  length  of  the  series  studied  by  Hurst.  However, 
by  returning  to  equation  4.27,  we  note  that  taking  a  and  b  so 

i, 

that  b  *  a  does  not  change  the  asymptotic  growth  rate.  We  now  note 
from  Eq,  4.24  that  in  this  case, 

mR  (a,b;t)  =  mR(l ,1 ;at)  . 

Hence  by  choosing  a  sufficiently  small,  we  equivalently  perform 
a  time  expansion,  so  that  the  approach  to  the  asymptotic  value  can  be 
delayed  as  long  as  desired.  Now,  if  our  model  is  to  display  the  Hurst 
phenomenon  up  to,  say,  t  =  2000  (years),  we  would  need  to  take 
a  3  .001.  However,  such  a  value  of  a  implies  that  the  correlation 
between  successive  yearly  increments  is,  from  Eq.  (3.3), 
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CorrC&j  (t)  ,A!  (t+1  )]  =  -j-1 - =  -99865 

21. 002-1+e  -I 

This  value  is  much  too  high  to  be  realistic.  However,  if  we  take 
a  =  .1  ,  then  we  obtain 

Corr[Aj (t) ,t\ (t+1 )]  =  .87720. 

Thus  a  stream  flow  with  a  yearly  correlation  of  this  value  and  following 
this  model  would  display  Hurst  behavior  out  to  approximately  20  years. 


rison  of  growth  rate  of  expected  range  for  Markov  model  to  average  growth  rate 
irst  data  for  *1  <  t  <  2  . 


CHAPTER  V 


INVARIANCE  METHODS  IN  THE  ANALYSIS  OF  SEMI-INFINITE  AND  FINITE  DAMS 
5.1  General  Description 

In  this  chapter  we  will  develop  the  application  of  what  we  call 
in  general  “invariance  methods"  to  the  analysis  of  storage  theory 
problems.  These  methods  are  well-known  in  the  fields  of  astrophysics 
and  particle  physics,  where  the  power  of  these  techniques  have  been 
recognized  and  well  utilized.  Originally  formulated  by  V.  A. 

Ambarzumian  [3]  in  1944,  the  techniques  were  formally  developed  by  S. 
Chandrasekhar  [16]  in  a  series  of  papers  on  radiative  transfer  through 
semi-infinite  stellar  atmospheres  which  appared  in  the  Astrophysics 
Journal  from  1944  to  1947.  Subsequently  Bellman  and  Kalaba  in  1956  [5] 
refined  the  technique  by  applying  it  to  finite  regions  in  the  field  of 
particle  transport,  in  which  they  called  the  technique  'invariant 
imbedding.'  A  good  definition  of  invariance  methods  is  given  by  M. 

Scott  [40]  who  states  that  the  method  involves: 

generating  a  family  of  problems  by  means  of  a  single 
parameter,  where  the  basic  properties  of  the  system  remain 
invariant  under  the  generation  of  the  family.  The  family 
then  provides  a  means  of  advancing  from  one  member,  some¬ 
times  degenerate,  to  the  solution  of  the  original  problem. 

After  providing  a  brief  description  of  Chandrasekhar's  and  Bellman's 
application,  we  will  establish  what  will  be  a  natural  application  of 
these  methods  to  storage  theory.  In  particular  we  will  take  advantage 
of  the  power  of  the  invariance  methods  to  derive  the  Laplace  transform 
of  the  'wet  period,'  and  from  this  obtain  moments  and  the  limiting 
probability  of  emptiness.  This  will  lead  us  to  a  discussion  of  neces¬ 
sary  and  sufficient  conditions  for  the  recurrence  of  the  contents 
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process.  Although  some  of  the  results  that  we  present  in  this  section 
are  not  new,  we  will  be  able  to  derive  them  in  a  direct  and  simple 
fashion,  avoiding  the  Kolmogorov  differential  equations  which  have 
previously  been  used  to  derive  these  results.  In  fact,  in  the  case  of 
the  semi-infinite  dams  our  equations  are  algebraic  in  nature. 

5.1.1  Chandrasekhar's  Principles  of  Invariance 

In  his  application  to  astrophysics,  Chandrasekhar  was 
concerned  with  the  transfer  of  radiation  through  a  stellar  atmosphere. 

In  studying  certain  aspects  of  the  radiation,  such  as  the  intensity,  he 
considered  the  atmosphere  to  be  stratified  in  parallel  planes  in  which 
all  of  the  pertinent  physical  properties  are  invariant  over  a  plane. 

When  considering  a  semi-infinite  atmospheric  region  then,  he  utilized 
the  fact  that  certain  physical  laws  governing  radiative  transfer  must 
be  invariant  to  the  addition  or  subtraction  of  layers  of  arbitrary 
optical  thickness  to  or  from  the  atmosphere.  Through  this  'principle 
of  invariance'  as  he  called  it,  he  was  able  to  increment  the  arguments 
in  the  equation  governing  the  physical  laws  of  radiative  transfer 
while  still  maintaining  the  basic  equations.  These  perturbations  of  the 
equations  then  led  him  to  direct  solutions  for  the  physical  quantities 
of  interest. 

In  studying  the  time  to  emptiness  in  a  semi-infinite  dam,  the 
essential  similarities  of  the  two  problems  become  clear.  In  the  case  of 
the  dam,  the  strong  Markovian  character  of  the  bivariate  process  (Ct,Xt) 
makes  it  possible  to  consider  the  process  as  regenerative  at  the 
stopping  times  corresponding  to  first  entrance  into  an  upper  semi¬ 
infinite  region  after  traversing  an  incremental  slab  of  arbitrary 
thickness.  Thus  random  variables  such  as  first  passage  times  become 


'invariant1  to  the  addition  or  subtraction  of  these  finite  incremental 
slabs.  Seen  in  this  context,  one  notes  the  similarity  between 
Chandrasekhar's  Principles  of  Invariance  in  stellar  atmospheres  and  the 
regenerative  property  of  a  Markov  process.  In  fact,  this  entire  dis¬ 
cussion  could  be  phrased  in  terms  of  regeneration.  The  essential 
difference,  however,  is  that  in  the  classical  analysis  of  Markovian 
structures  the  principle  of  regeneration  is  used  to  establish 
Kolmogorov-type  equations  which  are  generally  partial  differential 
equations.  These  equations  normally  provide  much  more  information 
about  the  process  than  is  required  for  the  study  of  particular  aspects 
which  become  essentially  boundary  conditions  for  the  Kolmogorov  equa¬ 
tions.  For  example,  suppose  that  we  are  interested  in  the  wet  period 
for  a  semi-infinite  dam:  that  is,  the  time  elapsed  from  an  initial 
exit  from  zero  to  a  first  return  to  zero.  In  the  two-state  case,  if 
we  let 

T  =  i nf { t  >  0  :  C^.  =  0) 

then  we  may  be  interested  in  finding  the  Laplace  transform  of  T  , 

*s  =  Ete~ST|Xo  =  ui  >  0 

Now,  by  setting  up  and  solving  the  Kolmogorov  backward 

equations,  Brockwell  was  able  to  determine  the  Laplace  transform  of  the 

time  to  first  emptiness  given  an  arbitrary  initial  level  x  .  Thus  if 

$  (x)  represents  this  transform  in  the  general  case,  then 

- 

♦ls(x)  ■  E[e'ST  |xo  ■  »1 
where  T  *  inf{t:t  >  0  ,  Ct  *  -x} 
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However,  the  solution  depends  on  <j> ^ s ( 0 )  ,  which  in  Brockwell's 
formulation  was  obtained  from  a  rather  complicated  set  of  boundary 
conditions.  Invariant  imbedding  leads,  as  we  shall  see,  to  a  direct 
determination  of  $  as  the  solution  of  an  initial  value  problem  when 
the  dam  is  finite  and  of  an  algebraic  equation  when  the  dam  is  infinite. 

5.1.2  Bellman's  Invariant  Imbedding 

As  noted  earlier,  8ellman  refined  Chandrasekhar's  Principles 
of  Invariance  and  applied  them  to  finite  regions  in  particle  transport 
theory.  Two  excellent  comprehensive  studies  on  these  applications  are 
those  by  Bellman,  Kalaba,  and  Prestrud  [6],  in  which  they  apply  the 
technique  to  radiative  transfer  in  slabs  of  finite  thickness,  and 
Bellman,  Kalaba,  Prestrud,  and  Kagiwada  [7],  in  which  they  apply  it  to 
time  dependent  transport  processes. 

In  these  studies.  Bellman  called  this  technique  "invariant 
imbedding,"  a  name  by  which  it  has  subsequently  become  known.  The 
application  to  time  dependent  transport  processes  is  particularly 
relevant  to  our  study  of  the  finite  dam.  Because  the  analogy  to  the 
finite  dam  problem  is  not  obvious,  we  will  discuss  briefly  the  former 
process,  and  then  draw  the  analogy. 

The  physical  model  in  the  transport  process  is  of  a  one-dimensional 
rod  of  fixed  length  z  which  is  capable  of  transporting  particles  such 
as  neutrons.  The  particles  are  allowed  to  travel  to  the  right  or  to 
the  left  and  can  interact  only  with  the  fixed  constituents  of  the  rod. 
When  a  particle  interacts  with  the  rod,  the  old  particle  disappears 
and  two  new  ones  appear,  one  travelling  to  the  right  and  one  travelling 
to  the  left.  All  the  particles  travel  with  the  same  speed  and  their 
other  physical  properties  are  such  that  the  particles  are  distinguishable 
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only  by  their  direction.  Suppose  that  we  inject  one  particle  at  one 
end  of  the  rod,  and  are  interested  in  the  total  number  of  particles 
emanating  from  the  other  end.  One  approach  to  this  problem  would  be 
to  derive  an  equation  for  u(x)  =  number  of  particles  emanating  from 
the  end  of  the  rod  if  we  start  with  one  particle  at  x  ,  for  0  <_  x  <  i  . 
The  desired  quantity  would  then  be  given  by  u(o)  .  However,  a  more 
efficient  approach  would  be  to  develop,  if  possible,  an  equation  for 
u(o)  directly.  This  is  precisely  the  approach  used  in  invariant 
imbedding.  The  invariant  imbedding  technique  is  to  embed  the  original 
rod  of  length  i  into  a  rod  of  length  l  +  Si  ,  and  develop  a  dif¬ 
ferential  equation  for  u(t)  .  Whereas  a  differential  equation  for 
u(x),  0  <_  x  <  i  would  lead  to  a  two-point  boundary  value  problem  with 
boundaries  at  x  =  0  and  x  -  i,  the  differential  equation  for  u(o) 
is  an  initial  value  problem  with  initial  value  l  -  0  .  This  initial 
value  is  easily  determined  because  of  its  degeneracy. 

We  now  begin  to  see  the  analogy  to  the  storage  problem,  which  we 
will  list  in  the  form  of  a  table. 


Particle  Transport 


Storage  Problem 


1.  Particles  travel  only  to  right 
or  left  with  fixed  velocity. 

2.  Do  not  need  to  know  entire 
function  u(x),  0<x<i,  but  only 
u(o). 

3.  Differential  equations  for 
u(x)  lead  to  a  two-point 
boundary  value  problem,  where¬ 
as  differential  equations 

for  u(o)  give  an  initial 
value  problem. _ 


1.  Dam  contents  only  increase  or 
decrease  at  fixed  rates. 

2.  Do  not  need  to  know  (x)  for 
all  0  <_  x  <_  a  ,  but  onTy  i>s(o) 

3.  Differential  equations  for 

<t> -  (x)  lead  to  a  two-point 
boundary  value  problem,  where¬ 
as  differential  equations  for 
<P  (o)  gives  an  initial  value 
problem, _ 


Table  5.1  Analog  between  Particle  Transport  and  Storage  Problem. 


68 


5.2  Applications  of  Invariance  to  First  Passage  Times 

We  will  now  use  invariance  techniques  to  study  the  distribution 
of  first  passage  times.  In  a  later  section  we  will  discuss  necessary 
and  sufficient  conditions  for  these  times  to  be  finite. 

5.2.1  First  Passage  Times  for  the  Semi-Infinite  Dam 
5. 2. 1.1  General  Case 

We  define  the  following  functions  related  to  first 


passage  times: 

(a)  The  Laplace  transform  of  the  first  passage  time  to  level  y  : 


Let  t  =  inf{u  :  Cu  =  y}  , 


and  define 


Tij  ’  E 


-St 


'(*  •  ,.>!xo  ’  “i 

T  J 


(5.1) 

(5.2) 


for  u-  >  0  ,  Ui  >  0  ,  y  >  0  and 

“  E  [e"S',<xT  *  Uj>lX°  *  “']  (5-3) 

for  u1-  <  0  ,  u,.  <  0  ,  y  <  0  . 

If  there  are  p  positive  states  and  m-p  negative 
states,  let  T+  and  T~  represent  the  corresponding  matrices 
with  entries  T+.  and  TT .  ,  of  dimension  pxp  and 

•  J  1  J 

(m-p)  x  (m-p). 

(b)  The  Laplace  transform  of  the  first  return  time  to  level  zero: 


For 


u,j  >0  ,uj 


(5.4) 


where 
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T  =  inf{t  >  0  :  Ct  =  0} 


Let  R+"  be  the  px(m-p)  matrix  C R"t T ]  . 

For  Ui  <  0  ,  uj  >  0,  the  right  side  of  (5.4)  will  be  denoted 

by  RTj  and  the  (m-p)xp  matrix  [R”p  by  R~+  . 

(c)  The  Laplace  transform  of  the  first  return  time  to  level  zero 

with  no  prior  passage  through  a  given  level. 

For  Ui  >  0  ,  u-  <  0  ,  and  y  >  0  ,  let 

+  r  -sT  i 

Rij(y)  =  EL6  !{xt  =  u-l^CU^y,  0<  t<  T  >iXo  =  ui 

J  "  J(5.6) 

For  u.  <  0  ,  u  >  0  ,  let 

Rij(y)  =  E  e  hx-  =  vA  !{C(t)^-y  0<t<T  }lXo  =  ui 

1  J  ~  (5.7) 

Let  R+(y)  be  the  px(m-p)  matrix  [Rt-(y)]  and  let 

R  (y)  be  the  (m-p)xp  matrix  [R^j(y)]  > 

The  Laplace  transform  of  the  first  return  tine  of  the  content  of 
a  topless  dam  to  zero  given  that  XQ  =  is 

l  <  j  • 

jiu-O  J 

where  is  defined  by  (5.4).  We  now  use  an  imbedding  argument  to 

determine  R+~  .  The  time  at  which  a  sample  path,  starting  at  the  zero 
level,  crosses  into  the  upper  region  (y,«)  is  a  stopping  time,  and 
the  strong  Markov  property  of  (Xt,Ct)  may  be  applied  at 
these  times.  Using  this  idea  to  obtain  the  desired  first  return  time 
to  zero,  we  decompose  R  according  to  the  number  of  reflections  of  the 
sample  path  in  the  slab  (0,y)  to  obtain  the  relation 
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R+-  =  R+(y)  +  T+(y)R +"T'(y) 

+  T+(y)R+'R  "(y)R+"T'(y) 

+  T+(y )R+"R  “(y)R+‘R  '(y)R+"T‘(y) 

+  •  .  • 

=  R+"(y)  +  T+(y)  R+'[l  -  R  ”(y)R+']“1T'(y ) 


We  now  consider  the  limiting  behavior  of  the  matrices  in  (5. 
y  .  First,  note  that 

T|i(y)  =  ^l-xi  e  yi  +  o(y)  as  y-0 
and,  expanding  the  exponential, 

Tii(y)  =  1  ■  -f:  (xi+s)  +  o(y) 

Similarly, 

T|j(y)  =  Vij  ^  e  Wi  +  0(y)  ’  ii4j 


-  Vy^T 


+  o(y) 


Therefore, 

T+(y)  *  I  -  y  [d_1(sI-Q)]  ++  +  o(y) 

where,  for  an  mxm  matrix  M  we  define  [M]++»  [M]+_  ,  [M]_. 
[M] _  to  be  the  partitions  defined  by 


M 


P 

(m-p) 


(5.8) 

(5.9) 

for  small 

(5.10) 

(5.11) 

(5.12) 

(5.13) 

(5.14) 

,  and 

(5.15) 


P  (m-p) 
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Through  parallel  reasoning, 

T'(y)  =  I  +  y  [d~1  (sI-Q)]  __  +  o(y) 
For  the  R  functions,  it  can  be  shown  that 


(5.16) 


and 


R+(y )  =  y[D_1Q]+_  +  o(y) 
R"(y)  =  -y[D_1Q]_+  +  o(y) 


(5.17) 


(5.18) 


Substituting  (5.14),  (5.16),  (5.17),  and  (5.18)  into  (5.9)  we  obtain 
R+_=  y(D_1Q)+j-  ^i-yfo"1 (si-Q)j  ++j 
+  R+*Ll+y(0"  Q).+RfT2(i+y(o‘!(si-Q)  __j+  0 


(y) 


Since 


then 


f^y(0_1Q)_+R+' 


"l=  I-y(D_10)_+  R+*  +  o (y), 


(5.19) 


(5.20) 


(5.21) 


R+'  =  y(D_1q)+_  +  R  ‘  -  yjj)”1 (sI-Q)J  ++R  ' 

-  R+*y(0":Q)_+R+'  +  R+'y[D"l(sI-Q)]__  +  o(y) 

from  which  we  obtain,  by  taking  the  limit  as  y-*0  , 

r+'(D’1Q)_+  R+-  +[d_1 (sI-Q)J++  R+-  -  R+"[d_1 (sI-Q)]__ 


-  (D"lQ)+_  -  0  (5.22) 

5. 2. 1.2  Two-state  Case  -  Laplace  Transform  of  Wet  Period 
For  the  two-state  case  we  remarked  earlier  that  the 
first  passage  time  corresponds  to  the  wet  period  in  the  dam.  In  this 


72 


case  we  set  4>s  =  R12  in  (5.4)  and,  making  the  appropriate  substitutions, 
equation  (5.22)  reduces  to 


PM"1  A2  T  _1  ,  _1  r-t  _1  “I' 

2  <fr  -  Lf  --i-  - 


,s  .PU2  -Xu^  -s(y"  -u"  )]  *s-X  =  0 
from  which  we  obtain  the  two  solutions 


(5.23) 


4>s  =  M2  (2p ) 


”X|a  +  (A  +4Xpyj1u21 )** j 


(6.24) 


where 


.  ,  -1  -1.  ,  -1  -1 

A  3  S(U2  "Ml  )  +  PM2  "Aux 


We  must  now  select  the  correct  solution.  Since 
A2  -»•  «  as  s  -*•  <=  ,  then  ()>s  N 

is  unbounded  for  the  positive  solution,  and  so  we  conclude  that 
•  (2p )  1  A  -  (A  +4Xpui 1 


-  M2 ' 


(5.25) 


We  now  investigate  conditions  under  which  the  first  return  to  zero, 
T  ,  is  finite.  First  note  that 

-1 


and 


$0  =  M2  (2p  )  |^PU2 1  "  Xui1  -  j  p y 2 1  + 

PM2  +  Xyi  3  m  ^ 


(5.26) 


^M 1 M  2  /  > 

where  m  =  EXt  when  Xt  has  its  stationary  distribution.  We  call  m 
the  drift  of  the  process.  Therefore, 


1  if  n<0 
x  uz 

-  -  —  if  m  >  0 
P  Ml 


(5.27) 
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Since  4»  =  P[T  <  «]  ,  then  from  (5.27)  it  follows  that  for  the 

topless  dam,  a  non-positive  drift  implies  that  P[T  <  °°]  =  1  ,  whereas 

Au2 

a  positive  drift  implies  that  PL T  <  »]  - - 

i-  * 

5. 2. 1.3  Distribution  of  Wet  Period  for  the  Symmetric  Case 
If  we  set  m  =  b,  u2  =  -b»  x  =  o  =  a  ,  then  the 
Laplace  transform  equation  becomes 


£  t2  - 
b  vs 


2a  .  2s 

T+  "F 


>s  +F  =  0 


or 


K 


2  (1  +  +  1  =  0 


(5.28) 


From  this  equation  we  see  the  rather  surprising  result  that  the  time  to 
first  emptiness  does  not  depend  on  the  states  on  which  the  Markov  chain 
is  defined. 

Solving  the  quadratic  equation,  we  obtain 


=  1  +  -  -  (  i2-  +  2s 


(5.29) 


which  can  be  written  in  terms  of  1/s  as 


>s  •  f 


+  i 


Expanding  (l  +  1°  Its  binomial  expansion,  we  obtain 


>-  ?(*)(¥)  ‘1  ♦! 


(5.30) 


(5.31) 


Performing  a  term-by-term  inversion  of  (5.31),  we  obtain  ti.e 
density  of  T  ,  say  fy(x),  x  _>0  as 
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For  the  distribution  function,  say 

x 

<Ux)  =  P[T  <  x]  =  /  fT(u)du 

o 

wc  can  integrate  the  power  series  term  by  term  to  obtain  the  series 


5. 2. 1.4  Time  Between  Overflows  for  the  Semi -Infinite 
Bottomless  Dam 

We  now  consider  a  semi-infinite  dam  with  a  top  but 
no  bottom.  This  may  in  fact  be  a  more  realistic  model  than  the  topless 
dam,  since  actual  dams  are  constructed  and  operated  in  such  a  way  that 
the  probability  of  overflow  is  higher  than  the  probability  of  emptiness. 

The  analysis  of  the  bottomless  dam  is  directly  analogous  to  the 
analysis  of  the  topless  dam.  Therefore,  we  will  be  brief  in  the 
following  derivations.  In  the  bottomless  dam,  if  we  measure  the 
contents  relative  to  the  top  of  the  dam,  which  we  can  set  equal  to  zero 
without  loss  of  generality,  then  for  the  contents  process  we  have 
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B  -t 

=  j  X*(u)du 

o 


(a. 35} 


where 


X*(u) 


3  if  C®  =  0  and  X 
li  u 


0 


(X  otherwise. 

D 

For  the  bottomless  dam,  =  0  represents  an  overflow  condition. 
Let  0s  be  the  Laplace  transform  of  the  time  to  first  overflow,  given 
a  starting  condition,  XQ  =  p2  •  Then  by  considerations  of  symmetry  we 
can  conclude  that  satisfies  equation  (5.25)  with  the  parameters 
\,c,  and  u],u2  interchanged.  Therefore  }  satisfies 


,  -i  2 
Xui  <frs 


Ay/-  OU2  1  ~  S(u2  1  "U1  1  ) 


-  pM2  ‘  =  0 


(5.36) 


from  which  it  follows  that 


K  =  l  i  (2x )' 


A-(A2+4\oui1  U21)*5 


(5.37) 


where 


-l  _i .  _i  -l 

ft  =  S(m  -U2  j  +  -  P  U  2 


The  explicit  inversion  for  the  symmetric  case  is  identical,  since 
in  that  case  the  two  exchanged  parameters  are  equal.  Therefore,  the 
distribution  of  the  time  between  overflows  in  the  bottomless  dam  is 
equal  to  the  distribution  of  the  wet  period  in  the  topless  dam,  as 
expected. 

By  similar  reasoning  to  that  used  to  establish  (5.27),  we  conclude 
that  for  the  bottomless  dam,  a  non-negative  drift  implies  that 
PC T  "  “>]  =  1  ,  whereas  a  negative  drift  implies  Chat 

py  l 

P[T  <  -]  •  -  - — 

AU2 


5.2.2  First  Passage  Times  for  the  Finite  Dam 


5.2.2. 1  General  Case 


We  now  turn  our  attention  to  the  finite  dam  with 


an  upper  boundary  at  a  .  In  the  same  vein  as  (5.4),  we  define 


+  \  'ST 

Va)'EL;  1 


iXT  =  yj]'X°  '  Ui 


(5.38) 


p 

with  { Ct }  replaced  by  {Ct}  ,  the  content  process  for  the  finite  dam 

where  T  is  given  as  in  (5.5).  Let  R(a)  be  the  px(m-p)  matrix 

R+'(a)  . 

L  ij 

Using  the  same  matrices  defined  for  the  semi-infinite  dam,  and 
using  the  same  argument  of  decomposing  R+  (a)  according  to  the  number 
of  reflections  of  the  sample  path  in  the  slab  (0,y)  ,  0  <  y  <  a  ,  we 

can  write  the  relations 

cr 

R+'(a)  -  R+(y)  +  E  T+(y)R+’(a-y)  [V(y)R+'(a-y)]  V(y) 
i=U  L  J 

=  R+(y)  +  T+(y)R+_(a-y)  [l-R“(y)R+“(a-y)] -1T"(y) 

(5.39) 

Using  the  relations  (5.17)  -  (5.20),  we  tan  write 
R*"(a)  =  y(D-1Q)+_  +  (l-y[D_1 (sI-Q)] 

•  R+’(a-y)  [l-y(D_1Q)_+R+”(a-y) J  ( I+y [o- 1 ( s I -Q )]  )  +  o(y) 

(5.40) 

so  that 

R+'Ial  -  Rl'(a-/1  =  .  R+‘(a-y)(D_1Q)_+R+"(a-y) 

-[o-1 (sI-Q)]  ++R+"(a-y)  +  R+‘(a-y)[D"1 (sI-Q)] 
M0_1Q)+.  +  (5.4!) 


(5.41) 
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Taking  the  limit  as  y-*0  we  obtain 

d-d^  +  R+'(a)(D_1Q)_+R+'(a)  +  [D-1  (sI-Q)]++R+‘(a) 

-  R+"(a)[D_1 (sI-Q)] _  -  (D'1Q)+_  =  0  (5.*.  ) 

The  initial  condition  is 
R  (0)  *  Ee  Srij 

★ 

where  t-jj  is  the  first  passage  time  from  state  i  to  state  j  in 
the  Markov  chain.  This  variable  is  discussed  in  Chung  [17]. 

5. 2. 2. 2  Two-state  Case  -  Laplace  Transform  of  Wet  Period. 

In  the  two-state  case,  letting  $(a,s)  =  , 

equation  (5.42)  becomes 

+  ou2l«2(a,s)  +  [uj1  (a+s)  -  u^1  (p+s)]  $>{a,s)-\u^13  0 

(5.43) 

The  initial  condition  in  this  case  comes  from  the  holding  time  in  state 
ui  ,  which  is  exponential  with  parameter  x  .  Hence  we  have  the 
initial  condition 

<t>(0,s)  =  x(x+s)'i 

We  will  carry  out  the  complete  solution  for  this  case. 

Equation  (5.43)  is  a  Riccati  equation,  which  can  be  reduced  to  a 
second  order  linear  equation  by  letting 

4>(a,s)  =  y'(a,s)/bQy(a,s) 
where  bQ  *  pua1 

b-j  =  ui 1  (a+s)-^1  (p+s) 

=  ■  Am i 1 


(5.44) 
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The  resulting  2nc*  order  linear  equation,  as  can  be  easily  checked,  is 
y"(a,s)  +  b1y'(a  ,s )  +  b0b2y(a,s)  =  0  (5.45) 

To  solve  (5.45)  we  first  obtain  the  characteristic  roots: 


uz  (p+sj-m1  (x+s) 

0  a  - r -  + 


U21 (p+sj-ui1 (X+S) 


+  pXui ^  U2 


(5.46) 


Note  that  the  discriminant  is  equal  to 


U22(p+s)2  U1  (X+s)2 


-  Ml  1p2ls(p+X)-Ul 


1  -1 


(5.47) 


4  4 

which  is  always  positive  because  U2 <  J  •  Therefore  the  characteristic 
roots  are  real,  and  the  solution  for  y(x,s)  is: 


6xa  02a 
y(a,s)  =  Cxe  +  C2e 


Substituting  (5.48)  into  (5.44)  we  obtain 

exa  020 

C 1 9 i S  +  C2©2® 

i>(a,s)  = 


-1 

PU2 

0  1  d  9 

Cxe  +  C2e 

Now,  letting  C  =  C2/C x  ,  we  obtain 


$(a,s) 


exa  e2a 

9xe  +  C02e 


J  0ia 

M  3  + 


@2® 


(5.48) 


(5.49) 


(5.50) 


pu2  1  +  Ce 

which  involves  only  one  arbitrary  constant.  To  evaluate  this  constant 
we  use  the  initial  condition  to  obtain: 
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4>  ( 0 ,  s ) 


A 

A  +  S 


8j  +  C  0i 


PU2"1 0+c) 


(5.51) 


Hence 


(A+sJQj  -  0V2  1 ^ 


C  =  - 


-1  , 


(A+S  ) 0 2  '  PU2  x 


(5.52) 


So  we  have  as  a  final  result  that 


<t>(a,s)  = 


0  j  cl  02^ 

@ie  +  C82e 


PU2' 


9 1  a  0  2a 
e  +  Ce 


(5.53) 


where  C  is  given  by  (5.52)  and  e1,e2  are  given  by  (5.46). 


5.3  Moments  of  Time  to  First  Emptiness  for  the  Two-State  Case 

Since  we  have  derived  the  Laplace  Transom  of  the  time  to  first 
emptiness,  we  can  use  these  equations  to  obtain  the  moments  without 
having  to  actually  solve  for  the  transform.  Although  we  could  in 
theory  write  down  the  moment  equations  for  the  general  case,  the 
resulting  equations  are  not  easily  solved,  and  so  we  restrict  ourselves 
to  looking  at  the  first  two  moments  in  the  two  state  case,  for  which 
the  equations  are 


ET 


d_ 

ds 


s=0 


ET 


2  -  Ji 


ds2 


ps 1 s=0 


so  that 


VjrT-^,sL-o  ‘  I.5T  ♦,] 


1  s=0 
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For  the  semi-infinite  dam  we  differentiate  equation  (5.23)  to 
obtain 

2d  v>2  1 ( S  )<{>  (s)  -  [pu21_  AUi*"  S  (u/-  M21)]^c  (5  .D') 


,  -1  -1 .  „ 
+  (mi  -  m2  )$s  =  o 


And,  upon  evaluating  this  expression  at  s  *  0  we  obtain 


-2pu2  ET  +  (pw 2  l~  A.U1  )ET  +  Ml  -U2  =  0 


from  which  we  obtain 

M2  -  Ml 


ET  - 


U10+U2* 

which  can  be  written  in  terms  of  the  drift  parameter  as 


/  1  \  I  M2-M1  , 

ET  "  m-° 


(5.55) 


(5.56) 


Differentiating  equation  (5.23)  again,  we  obtain 


2py2J£>(s)i|>"(s)  +  r(s)2]  -  ^pm21-Xmi1-s(mi1"M21  )]  <t>' 


+  2[mi1-m21]((i's  =  0 


(5.57) 


Evaluating  again  at  s  =  0  we  obtain 


2py;l[ET2+E2T]  -  [puil-Xur1]ET2-2(yrI-ui1)ET  =  0 
From  which  it  follows  that 

Z(i‘1-V21)eT  -  Zpu^E2! 


-l  .  -l 


-l  -i ' 


ET2 


-l  -l 

PU2  “  Au2 


(5.58) 


Therefore 


Var  T  = 


2(m  U2  )Ef  -  (3pp2  -Ayj  )E“T 


p  +  A 


(5.59) 


Similar  calculations  for  the  finite  dam  produce 


ET 


+  PU2  ) 


-1 

Hi 


-Uu^W  )a ' 

)e 


(5.60) 


5.3.1  Conditions  for  Recurrence 

In  analogy  with  the  theory  of  Markov  chains,  we  call  the 
contents  process  Ct  recurrent  or  transient  depending  on  whether 
P  -  1  or  p  <  1  where  p  -  P[T  <  »>]  .  In  the  recurrent  case  we 
call  Ct  positive  recurrent  or  null  recurrent  depending  on  whetner 
FT  <  ®  or  Ej  *  »  ,  We  will  relate  conditions  for  recurrence  to 
the  drift  parameter  of  the  process,  m  . 

For  the  topless  dam,  it  follows  from  equation  (5.27)  that  m  <_  0 
implies  recurrence,  while  m  >  0  implies  transience.  From  equation 
(5.56)  it  follows  that  m=0  implies  null  recurrence  while  m  <  0 
implies  positive  recurrence. 

By  symmetry,  for  the  bottomless  dam  it  follows  that  m  >  0  implies 
positive  recurrence,  m=0  implies  null  recurrence,  and  m  <  0  implies 
transience. 


5 . 4  Limiting  Probability  of  Emptiness  for  the  Two-State  Case 
Consider  the  time-dependent  probability  of  emptiness 

Pe(t)  =  P(Ct  =  0)  . 

If  the  bivariate  process  has  a  limiting  distribution,  then  we  consider 


4 
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the  following  1 imit 


P 

e 


Tim  Pg(t) 


We  call  this  the  limiting  probability  of  emptiness,  and  it  can  be 
interpreted  as  the  large-time  percentage  of  time  that  the  dam  will  be 
found  empty. 

Brockwell  [.12]  has  actually  derived  the  limiting  distribution 
of  the  content  by  solving  the  Kolmogorov  equations.  However,  as  is  the 
case  with  the  time  to  first  emptiness,  these  equations  require  knowledge 
of  the  boundary  value  Pg  for  their  solution.  Although  Pg  can  be 
obtained  by  other  methods,  we  present  here  a  direct  method  of  evaluating 
Pg  by  using  the  previous  results  concerning  the  time  to  first  emptiness 
and  a  renewal  argument. 


Embedded  Regenerative  Process 


Assuming  that  X(0)  =  m  >  0  ,  define 


ti  =  i nf { t  >  0|C^  =  0} 
t2  =  i nf { t  >  tj |Ct  =  0} 


and,  continuing  in  this  fashion,  define  t.  ,  i  >_  1  .  The  x.'s  are 
the  first  return  times  to  zero  of  the  contents  process,  and  since  the 
bivariate  process  (Xt ,Ct )  is  strong  Markov  and  the  t's  are  stopping 
times,  then  they  are  regeneration  points  and  the  process  (tj.ti  +■  xz,-..) 
is  an  embedded  delayed  renewal  process  (see  Figure  5.1). 


Figure  5.1  Embedded  Renewal  Process. 


Now,  thinking  of  this  as  a  two-cycle  renewal  process  with 
cycles  T  and  X  ,  then  it  is  a  standard  result  of  the  renewal  theorem 
that 


P  S  - _ 

e  EX  ET 


5.4.2  Semi -Infinite  Topless  Dam 

For  the  semi-infinite  topless  dam  we  have 


and,  by  the  memoryless  property  of  the  exponential, 


(5.61) 


(5.62) 


X  -  exp  (? ) 
so  that  EX  =  n"1- 


(5.63) 


Substituting  (5.62)  and  (5.63)  into  (5.61)  we  obtain  after  some 


5.4.3  Semi-Infinite  Bottomless  Dam  -  Limiting  Probability  of  Over¬ 
flow 

Using  the  symmetry  between  the  topless  and  the  bottomless 
version,  we  can  write  for  the  time  to  first  overflow,  T 


and,  using  the  analogous  embedded  renewal  process,  we  have  for  the 
limiting  probability  of  overflow,  , 


p  =  _ il! _  =  JL 

X  +  Ui^+UlP 


(5.66) 


5.4.4  Finite  Dam 


In  this  case,  substituting  equations  (5.60)  and  (5.63)  into 
(5.61 )  we  obtain 


Pp  3  +  PU21] 


(X+p)(Xyi 1  +  PU21® 


-(xui1  +  pu^Mal 


(5.67) 


For  the  case  ui  =  1 ,  U2  3  -  1  we  obtain 


Pe  3  x(p-x) 


/  <p-A)a  u-1 

U+p)\pe  -  x 


(5.68) 


which  agress  with  the  previously  reported  result  [12]. 


CHAPTER  VI 


A  FUNCTIONAL  CENTRAL  LIMIT  THEOREM  FOR  THE 
UNRESTRICTED  CONTENTS  PROCESS  ON  TWO  STATES 

6.1  Formulation 

In  this  chapter  we  show  that  a  scaled  version  of  the  contents 
process  Ct,  0  <_  t  <  ®  ,  will  converge  weakly  to  the  Wiener  process  on 
D[0,*>).  A  short  description  of  the  topology  of  this  space  was  pro¬ 
vided  in  Chapter  II.  The  weak  convergence  results  of  this  section 
generalize  the  results  of  Fukushima  and  Hitsuda  [22]  and  Pinsky  [33], 
who  showed  convergence  of  the  marginal  distributions. 

To  establish  weak  convergence  for  Ct  we  appeal  to  a  result  given 
in  Chapter  II  that  if  Ti  are  the  times  of  jumps  of  a  Markov  chain 
and  (Xp }  are  the  succession  of  states  visited,  then  if  the  sequence 
{ }  is  deterministic,  the  holding  times  are  unconditionally  independ¬ 
ent,  i.e., 

n  -\  .  ,u  . 

P[T . -T .  >  u,,j  =  l,...,m]  =  it  e  J~  J 

J  J-i  J  j»i 

In  the  two-state  case  the  sequence  Xn  is  deterministic,  being 
given  by  (u i *U2»ui >U2»*  •  •  J*  If  we  fix  on  the  entrance  times  into  the 
ui  state  by  the  process,  we  see  that  we  can  consider  an  embedded 
renewal  process  with  cycle  lengths  (5^)  where 

50  =  0 

51  "  T2i  "  T2 ( i - 1 ) 
as  seen  in  Figure  6.1. 
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From  Figure  6.1  we  can  see  that  the  embedded  renewal  process  is 

really  an  alternating  renewal  perocess  with  two  subcycles  in  each 

interval.  Thus  the  first  cycle  has  length  £i  ,  and  is  composed  of 

the  two  subcylces  of  length  Ti  -  exp(A)  and  -  Ti  -  exp(p)  .  The 

second  cycle  has  length  £2  and  is  composed  of  the  two  subcycles  of 

length  T3  -  T2  -  exp(x)  and  T4  -  T3  -  exp(p),  and  so  forth.  The 

slopes  of  the  sample  path  over  odd  subcycles  is  u 1 ( >0 )  and  over  even 

subcycles  is  U2(<0)  .  We  use  the  notation  W*  to  represent  the  wait- 

t  h 

ing  time  random  variable  for  the  subcycle  of  the  i  n  cycle  during 
which  the  slope  of  the  sample  path  is  positive,  and  for  the  sub¬ 
cycle  during  which  it  is  negative.  Thus 


and 


“i ' 


u  =  T  -  T 
wi  2i  2 i - 1 


Note  that  W.  +  WT  =  T£i  -  T2(i_1}  =  ^  .  Let  (S^.,  be  the 
renewal  process,  so  that 
n 


Sn  =  Ci 


(6.1) 


Let  N(t)  be  the  renewal  counting  function  for  the  Sn  process;  thus 


N(t)  =  inf  {j |  [  i  t} 

i=l  1 


Now  consider 


nt 


C (nt )  =  /  x(s)ds 
0 


(6.2) 


(6.3) 


Then,  referring  to  Figure  6.1,  we  can  write  C(nt)  as 


38 


C(nt)  = 


f  sn(nt)  r 

J  x(u)du  -  J 

r\  S 


'5N(nt) 

x(u)du 


(6.4; 


nt 


where  the  second  term  represents  the  'overshoot'  of  the  process  to  tne 
first  renewal  past  nt  ,  and  so  must  be  subtracted  off  to  maintain 
equal ity. 

Since  the  sample  paths  of  X(t)  are  piecewise  linear  with  slopes 
Hi  and  uz  ,  we  can  re-write  C ( nt )  in  terms  of  the  holding  time 


variables  W.  as 

‘  N(nt) 


N(nt) 


C(nt) 


»  £  W*  Hi  +  £  WT 


**:  -  Rn(t) 


(6.5) 


where 


5N(nt) 


Rn(t)  =f  x  (u  )du 


nt 


Now  define 


+  Vi  U2 

wi  +  wi  U2  ‘  T  •  T 


(6.6) 


(6.7) 


Then,  by  Donsker's  theorem  (Billingsley  [9]), 


Ind. 

/n 


H?  u| 

w|  I  ~  + 

x2  P2  > 


on  D[0,») 


(6.8) 


Returning  to  the  original  process  (6.5),  we  can  write 
N(nt)  /  \ 

C  ( n  t )  L  (<u1+WTu 


+  lT+  t)  N(nt)  +  Rn(t) 


(6.9) 
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or, 

C(nt)  '  ZN(nt)  +  (t  *  T/  «’'l  *  V"  (6-,0» 

Setting  n  (t)  =  Mill  f  then  (nn(t)}“s^  are  increasing  processes 
in  t  ,  and  for  each  fixed  t  , 

P  t 

nn(t)  -*■  as  n-+®  . 

Therefore  it  follows  that 


nn( » )  =>  LE(c)]_1(*)  on  D[0,»),  i.e.. 


(6.11) 


Since  the  limit  is  degenerate,  then  the  following  joint  convergence 
holds  by  Theorem  4.4  in  Billingsley  [9]- 


Using  the  continuous  mapping 

i ^(x ,4> )  =  x  •  $ 
we  conclude  that 


LNin * ) J  — >  W(S’)  on  D[0,») 
Sn 


(6.12) 


(6.13) 
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We  now  show  that 


✓ff 


0  in  0(0,®)  .  From  6.6  we  have 


V'>  / 


5N(nt) 


x(u)du 


nt 


^(nt) 


/•M"  W 

Sup  x(u 

ntiulSN(n 

nt 


(u)  du 
N(nt) 


£  V(yj,-y2) 


=  V(uj ,-u2 )  y  (nt) 


where  y(nt)  is  the  forward  recurrence  time  evaluated  at  nt.  There¬ 
fore  it  suffices  to  show  that 


1^1  0  in  0[0,»)  . 

/if 


If  u,o2  are  the  mean  and  variance  of  the  interarrival  times  for 
a  renewal  process,  then  Iglehart  and  Whitt  [26.1]  have  shown  that 


.  N.Ln.* PJa  |  (w(.),  “W( • ))  .  (6.14) 

/na2/y3  J 


Now,  since  -ln. *■!■»> Ill 
n  u 


in  D[0,®)  ,  then  we  also  have  the  joint 


convergence 


uN(n*)-n- 

a/n 


Now  consider  the  mapping  ♦(x.y,#)  =  (x-<i>,y)  and  apply  the 
continuous  mapping  theorem  to  obtain 


t 


'■J*  ***»*■"? 
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l 


uN(n«)-  n- 
o  /n 


and,  again  applying  the  continuous  mapping  theorem  to 
4>(x»y)  =  x  +  y 

we  obtain 


(6.16) 


SN(n-)  ~  uN(n‘)  +  uN(n-)  -  n  • 
a  /FT 

°r  SWn-)  „ 

o/n 

so  that  indeed  ^ n =>0  on  D[0,»)  . 


6.2  Zero  Drift 

Ui  U2 

In  the  zero  drift  case  —  +  —  =  0  ,  so  that  from  6.10 

x  P 

it  follows  that 


C<”‘>  •  ZN(„t)  *  Rn(t> 
and,  using  Th.  4.1  in  Billingsley  [9]  , 


(6.17) 


”>  W(g»)  on  D[0,»)  . 

Hence,  if  AeS  ,  the  sigma  field  generated  by  the  open  sets  relative  to 
the  D[0,®)  metric. 


51ml  a 

/FT 


-  P  [w(B-)  eA  j 


Thus  we  have  the  large-n  approximation 


£ iSzl 
/n 


e  A 


=  P 


[w(bO  eaJ 


0  as  n-+» 
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6.3  Non-Zero  Drift 

For  non-zero  drift  Eq.  (6.10)  is 

C(nt)  ■  ZN(nt)  *  (t*  lr)N(nt)  *  Rn(t) 

so  that  if 

C*<„t>  ■  Slat)  .  Siail 

js  V  s/  js 

then 

C* (n - )— >  W(s-)  on  D[0,=°)  . 

Thus  we  see  that  in  this  case  the  centering  function  is  random. 


6.4  Asymptotic  Behavior  of  Range 


We  can  obtain  the  asymptotic  distribution  of  the  range  from  the 

corresponding  result  for  Brownian  motion  as  follows:  Consider  the 

continuous  mapping  from  D[0,«)  into  R  given  by 

f(c(nt))  =  ( flax  C(nu)  -  Min  C(nu)j  . 

\0<u<t  0<u<t  / 


Then  by  the  continuous  mapping  theorem,  for  R(nt),  the  range  of 

/n 

we  have  R 

R(nt)  ->  RB(et)  , 


the  range  of  Browian  motion. 

W.  Feller  [21]  derived  the  density  of  the  range  of  Brownian  motion 
as: 

f(t,r)  »  8 


where  <j> (x )  is  the  standard  normal  density  function.  Hence  the 
asymptotic  distribution  of  R(t)  has  a  density  given  by 


CHAPTER  VII 


SUMMARY  AND  CONCLUSIONS 


7.1  Summary 

In  this  investigation  we  considered  several  problems  connected  with 
Markovian  storage  models  in  continuous  time.  We  proceeded  by  con¬ 
sidering  three  different  variations  of  the  model.  These  were  the 
doubly-infinite  dam  without  top  or  bottom,  the  semi-infinite  dam  in 
both  its  topless  and  bottomless  versions,  and  the  finite  dam  with  both 
a  top  and  a  bottom.  In  each  case  the  process  which  we  examined  was 
the  bivariate  Markov  process  (Xt,  Ct)  in  which  Xt  was  a  Markov  chain 
with  finite  state  space  representing  the  net  input  rate,  and  Ct  was 
the  integral  of  Xt  ,  and  represented  the  dam  contents.  We  maintained 
as  much  generality  as  possible  throughout,  but  restricted  overselves  to 
particular  cases  whenever  necessary.  The  most  general  formulation  of 
the  model  was  when  the  Markov  chain  Xfc  was  defined  an  an  arbitrary 
n-dimensional  state  space  with  an  arbitrary  generator  matrix.  We 
often  restricted  the  discussion  to  the  two-state  case  in  which  the 
Markov  chain  was  defined  on  states  ^  >  0  and  uz  <  0  and  the  holding 
times  in  the  states  were  exponential  with  parameters  A  and  p  .  On 
occasion  we  restricted  ourselves  still  further  to  the  symmetric  case 
in  which  A  ■  p  and  v x  3  -  U2  • 

In  the  derivations,  we  did  not  follow  the  formal  definition-theorem 
style  of  presentation,  but  rather  chose  to  present  the  derivations  and 
results  in  a  more  natural  flow.  Since  with  this  style  of  presentation 
It  Is  not  as  easy  to  locate  individual  results,  we  now  present  the 

specific  results  derived  In  each  chapter,  as  an  aid  to  the  reader. 
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(iv) 


_2  2  - ( X+0 )s 

Cov[X^,X^.+  s]  =  pX(x+p)  ( y  i  “U2  )  e 

-(x+p)s 

Corr[Xt,Xt+s]  =  e 
Two-State  Case  -  Contents  process 


ECt  *  (p+X)  1(pwi+XM2)^ 


Var  3  2(p+x)  px(p 


Corr[Ah(t),Ah(t+h)]  3 


2[  -(x+p)t  1 

1-U2)  Le  +  (x+P)t-i  j 

.L 


■(x+p  )h 


2  (x+p)h-l+e 


■(X+p)h 


In  this  case  the  correlation  function  is  free  of  the  states. 


Chapter  IV 

In  this  chapter  we  studied  the  range  of  the  contents  process 
for  the  doubly  infnite  dam. 

(i )  General  Case 

We  derived  the  Laplace  transform  with  respect  to  time  of  a 
form  of  the  joint  distribution  of  the  maximum  and  minimum  variables. 

If  =  Sup{Cu,  0  <_  u  <_  T,  T  -  exp  (s)} 


and 


m$  »  Inf {Cu ,  0  <_  u  <_  T,  T  -  exp  (s)} 
•^(x.s)  3  1-P[M$  <_  a  -  x,  ms  >  -x|XQ  =  y.] 


1 


1 


then 

i(/(x,s)  3  expC-0*1  (Q-sI )x]  t|»(o,s) 


-  — - 
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where 

Va)  =  1,  ui  >  0 
yo)  -  1,  ui  <  0  . 


M$  is  given  by 
-ey 

e  ,  ui  >  0 

x  +  s  -  8uj  -ey 
- - -  e  ,  ^  <  0 


( i  i )  Two-State  Case 

The  distribution  of 

4  =■  yiui] = 


where 


s(*+o )+s2 
U1U2 


where 

A  -  P  *  s  +  *  *  s 
U2  Pi 


Hence  if  the  initial  rate  is  positive,  the  maximum  to  time  T  is 
exponentially  distributed.  If  the  initial  rate  is  negative,  the 
maximum  has  a  truncated  exponential  distribution  with  a  mass  at  0 

6lUl“S 

size  -j—  . 

If  R$  =  Ms  -  ms  ,  the  range  to  time  T  ,  we  found  that  when 
Xt  Is  stationary, 

ERS  *  (x+pT^er1  (p+x+s-ejm) 


+  9]'1[x+px'1(x+s+e1u2)] 


of 
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where  ex  =  e  above  and 
©!  =  9  -  A  . 


( i i i )  Symmetric  Case 

For  the  symmetric  case  in  which  yi  =  b,  v2  =  -t>,  and  \  =  p  =  a, 


we  found  for  the  range  to  time  t,  say  R  ,  that 


ERt-?  [F(-i-  ’■  -2at)  -  ’]  ■? 


-2a  t 


F  (  f-,  1,  2a t )  -  1 


where  F(a,B,x)  is  the  confluent  hypergeometric  series  defined  by 

V  (a)i  xJ 

F(a,3»x)  =  2-  7-f-  J T 

j=o  u'j 


where 


a(a+l  ). . . (a+j-1 )  for  j  =  l,2,... 


for  j=0 


Using  an  asymptotic  expansion  of  the  Laplace  transform  we  found 
the  asymptotic  behavior  of  the  range  as 

-3/2 


%(t)  - 


/E 


'z.  t  - 


2 

/a 


3  bt-j5_n[/2  b  t 

4/7  J72  t  -  0  3?  77 


r(-  1/2) 


so  that  in  particular  a  large-t  value  is 


mr(t)  ~  — 

r'  /a 


/* 


We  did  an  exact  calculation  of  ER^  for  0  <_  t  <_  10  in  the 

case  a  =  b  =  1  and  found  that  the  Hurst  phenomenon  was  evident  in 

the  approximate  range  0  <_  t  <_  2  .  If  a,  b  are  chosen  so  that 
!« 

b  =  a  ,  the  asymptotic  behavior  is  not  changed.  In  this  case,  choosing 
a  to  be  small  extends  the  Hurst  behavior  to  longer  periods  of  time, 
while  at  the  same  time  it  increases  the  correlation  between  increments. 


TOO 


In  order  to  explain  the  Hurst  phenomenon  with  this  model  for  Hurst's 
longest  series  of  about  2000  years  requires  a  correlation  between 
adjacent  yearly  increments  of  .99865,  which  is  unrealistically  high. 

For  the  shorter  time  period  of  20  years,  the  required  correlation  drops 
to  a  more  acceptable  value  of  .87720. 


Chapter  V 

In  this  chapter  we  introduced  the  principles  of  invariance 
developed  by  Chandrasekhar  and  Bellman  in  physical  models  to  the  study 
of  storage  models.  The  relevance  of  this  chapter  is  in  the  power  of 
these  techniques  to  provide  direct  solutions  to  first  passage  problems 
in  a  simplified  manner  over  classical  techniques. 

(i )  Semi-Infinite  Topless  Dam 
General  Case 

We  derived  Laplace  transform  equations  for  the  distribution 
of  the  first  return  to  zero,  given  an  exit  rate  >0  and  return  rate 

Two-State  Case 

In  this  case  the  first  return  time  to  zero  corresponds  to  the 
wet  period.  We  solved  the  general  system  of  Laplace  transform  equations 
explicitly  to  obtain  the  Laplace  transform  of  the  wet  period,  T  ,  as 

$s  =  U2 (2p )  l[A  +  (A2  -  4xpu i  Vz1)*5] 

where  A  =  sCu^-ui1)  +  PU2 1  ~  Xul1 


from  which  we  obtained 


P[T«»]  = 


i  if  m  ^  0 
phi 


Xy2 


if  m  =  0 


By  differentiating  $$  ,  we  obtained  (for  m  <_  0). 
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42"Ul 

0  A. 


Var  T  = 


2  { U 1  42  1  )E"*”  •  (3dM2  "Am  Me  T 

m(p+\) 


where  m  is  the  drift,  i.e.,  m  =  EX^  for  stationary. 

We  then  used  a  renewal  argument  to  obtain  the  limiting  probability 
of  emptiness,  Pg  as 

p  =  JL  (for  m  <  0) 
e  u2 

Synmetric  Case 

In  the  symmetric  case  we  were  able  to  invert  the  Laplace 
transform  to  obtain  the  density  of  the  wet  period  as 

fT(x)  =  |  F  (|,  3,  -2ax) 

where  F  is  again  the  confluent  hypergeometric  function. 

( i i )  Semi-Infinite  Bottomless  Pam 

In  this  case  we  derived  the  Laplace  transform  of  the  time 
between  overflows  as 

<t>s  =  ui(2a)  :[a  -  (A2+  4apui  Mz 
where  A  =  s(m  -  u2  1 )  +  Am  1  -  pu2 

The  density  of  the  time  between  overflows  is  the  same  as  the  density  of 
the  wet  period  for  the  topless  dam. 

We  derived  the  limiting  probability  of  overflow,  ,  as 


Pf  =  m/m  . 
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( i i i )  Conditions  for  Recurrence 

We  derived  necessary  and  sufficient  conditions  for  recurrence 
of  the  Ct  process  as  follows: 

Topless  Dam:  Ct  is  recurrent  iff  m  =  EXt  <  0  .  It  is  positive 
recurrent  if  m  <  0  and  null  recurrent  if  m  =  0  . 

Bottomless  Dam:  is  recurrent  iff  m  >_  0.  It  is  positive 

recurrent  if  m  >  0  and  null  recurrent  if  m  *  0  . 

(iv)  Finite  Dam 
General  Case 

We  derived  an  initial-value  differential  equation  for  the 
Laplace  transform  of  the  wet  period. 

Two-State  Case 

We  derived  explicitly  the  Laplace  transform  of  the  wet  period 

for  a  dam  of  height  a  as 

9  i d  9 

92e  +  C92e 

4>(a,s )  = 

where 

( X+S )9  j -py2  X 

C  - - -  ’ 

(X+S)02-PU2  \ 

9 1 ,2  =  A  +_  [A2+pXyi  lu2_1]  and 

uz  ^P+sJ-Ml  1  (A+S ) 

A  » - 2 - 

and  we  obtained  the  moment  of  T  as 
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-'  X-1  -1  -1.  -1  -l, 


ET  =  (\ui  +ou2  )  -U2  +  X  U2  (x+o)e 

and  the  limiting  probability  of  overflow  as 


•(Xu!  1+pu2  Ma 


pe  =  ^(p-a)  U+p  )  (pe 


(p-x)a 


Chapter  VI 


In  this  chapter  we  established  weak  convergence  in  0[0,»)  of 
the  unrestricted  contents  process  Ct  on  two  states.  The  basic  results 
which  we  obtained  involved  weak  convergence  to  Brownian  motion  of  the 
process  of  C(nt)  .  vie  established  that  for  the  zero  drift  case, 

-(n—*  -  w( b»  )  on  0[o,"),  where 


2  2 
P"U1  +  X^U2 


and  W(. )  represents  Wiener  measure. 

Using  this  we  were  able  to  conclude  that  the  asymptotic  distribu¬ 
tion  of  the  range  has  density 


7.3  Recommendations  for  Further  Study 

One  possible  line  of  further  research  would  be  the  extension  of 
the  results  derived  for  the  two-state  case  to  three  or  more  states. 
This  would  require  a  clever  choice  of  a  generator  matrix  which  would 
give  sufficient  generality  and  yet  have  enough  structure  to  make  the 
calculations  manageable.  One  possibility  which  we  would  recommend  for 
further  analyses  is  the  generator  matrix 
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Q  =  p(lir'-  I)  .  p  >  0 

which  has  the  property  that  tt  is  the  stationary  distribution  and  tlv 
autocorrelation  function  is 

„  ,  -ps 

Corr[Xt,Xt+s]  =  e  '  . 

Another  possible  line  of  research  is  the  extension  of  the  technique 
of  invariance  to  obtain  further  results  in  storage  theory.  These 
techniques  have  been  highly  developed  in  the  physical  fields.  What  we 
have  presented  here  is  but  the  briefest  introduction  to  their  use  in 
storage  theory.  Undoubtedly  many  other  results  must  be  obtainable 
through  sufficiently  clever  applications  of  these  principles. 
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